
==== Front
iScience
iScience
iScience
2589-0042
Elsevier

S2589-0042(24)01900-X
10.1016/j.isci.2024.110675
110675
Article
Deciphering sphingolipid biosynthesis dynamics in Arabidopsis thaliana cell cultures: Quantitative analysis amid data variability
Osinuga Abraham 1
González Solís Ariadna 23
Cahoon Rebecca E. 2
Alsiyabi Adil 14
Cahoon Edgar B. 2
Saha Rajib rsaha2@unl.edu
15∗
1 Department of Chemical and Biomolecular Engineering, University of Nebraska-Lincoln, Lincoln, NE 68588, USA
2 Department of Biochemistry and Center for Plant Science Innovation, University of Nebraska-Lincoln, Lincoln, NE 68588, USA
3 Center for Quantitative Cell Imaging, University of Wisconsin-Madison, Madison, WI 53706, USA
4 Natural and Medical Sciences Research Center, University of Nizwa, Nizwa, Oman
∗ Corresponding author rsaha2@unl.edu
5 Lead contact

05 8 2024
20 9 2024
05 8 2024
27 9 1106758 12 2023
22 4 2024
1 8 2024
© 2024 The Author(s)
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Summary

Sphingolipids are pivotal for plant development and stress responses. Growing interest has been directed toward fully comprehending the regulatory mechanisms of the sphingolipid pathway. We explore its de novo biosynthesis and homeostasis in Arabidopsis thaliana cell cultures, shedding light on fundamental metabolic mechanisms. Employing 15N isotope labeling and quantitative dynamic modeling approach, we obtained data with notable variations and developed a regularized and constraint-based dynamic metabolic flux analysis (r-DMFA) framework to predict metabolic shifts due to enzymatic changes. Our analysis revealed key enzymes such as sphingoid-base hydroxylase (SBH) and long-chain-base kinase (LCBK) to be critical for maintaining sphingolipid homeostasis. Disruptions in these enzymes were found to affect cellular viability and increase the potential for programmed cell death (PCD). Despite challenges posed by data variability, this work enhances our understanding of sphingolipid metabolism and demonstrates the utility of dynamic modeling in analyzing complex metabolic pathways.

Graphical abstract

Highlights

• Dynamic 15N labeling was used to estimate turnover fluxes of sphingolipids in Arabidopsis

• Developed regularized-DMFA framework effectively handles inter-sample variability

• Dynamic modeling reveals key enzymes crucial for sphingolipid homeostasis

Plant biology; Systems biology; Dynamic metabolic flux; Sphingolipid metabolism

Subject areas

Plant biology
Systems biology
Published: August 5, 2024
==== Body
pmcIntroduction

Sphingolipids represent a diverse class of lipids that mesh within the cellular structure of eukaryotic cells. Their significance is underscored by their dual role: serving as abundant component of cell membranes and bioactive molecules like ceramide (Cer) and sphingoid bases.1,2 With plants, typically, the significance of sphingolipids is undeniable. Constituting up to 40% of the plasma membrane, these lipids prominently feature in cellular endomembranes, such as the endoplasmic reticulum (ER), golgi, and tonoplast.3,4 Owing to their structural uniqueness, they have essential cellular functions. For instance, glycosylated variants such as glucosylceramides (GlcCers) and glycosylinositol phosphoceramides (GIPCs) not only bolster membrane function but also facilitate protein trafficking in the cell.2 Moreover, under environmental stress, such as microbial pathogenesis, an accumulation of specific sphingolipids triggers the onset of programmed cell death (PCD).5,6

With these molecules crucially situated in cell membranes and influencing cell signaling pathways,7 the sphingolipid biosynthesis (de novo), initiated in the endoplasmic reticulum (ER), is central to cellular function and homeostasis.7,8 The pathway specific to Arabidopsis thaliana (hereafter Arabidopsis), commences with serine palmitoyltransferase (SPT), catalyzing a foundational reaction involving palmitoyl-CoA and serine to produce 3-ketosphinganine, a rate-limiting phase in sphingolipids backbone, sphingoid, or long-chain base (LCB) synthesis.9 The ensuing enzymatic processes yield a variety of ceramides, further differentiated by ceramide synthases in Arabidopsis tailored to distinct acyl-CoA lengths.10,11,12 These ceramides, through glycosylation, form GlcCers, and further steps contribute to the synthesis of GIPCs, which are dominant glycosphingolipids in plant cells.4,13 The limiting and tight regulation, led by SPT,9 ensures balanced sphingolipid production while avoiding undue aggregation of PCD-inducing components.

In plant biology, the significance of sphingolipid biosynthesis is underscored by its regulatory responses to diverse stimuli, encompassing both biotic interactions and abiotic environmental factors. The interplay of genes, enzymes, and signaling molecules is fundamental to maintaining cellular homeostasis and facilitating adaptability to internal or external perturbations.14,15 Depending on the nature of perturbations—whether pathogenic or otherwise—cellular responses can manifest as transient dynamics or prolonged steady-state adjustments, reflecting the metabolic pathway’s adaptability across diverse physiological conditions. Notably, in plants, yeast, and human cells, the phosphorylation of sphingolipids and associated enzymes emerges as a pivotal, albeit not entirely deciphered, regulatory process.5,16,17,18

Studies suggest a connection between these phosphorylation events and signaling molecules, potentially triggering the production of reactive oxidative species, which are in turn implicated in PCD pathways.16,19,20 Because of this, enzymes such as long-chain-base kinase (LCBK) and sphingolipid-base hydroxylase (SBH) are noteworthy, as well as a putative sphingoid-phosphate phosphatase (LCBPP1) which is believed to dephosphorylate sphingoid molecules, restoring them to their long-chain-base precursor forms. Thus, these enzymes emerge as potential regulators in response to an external stimulus/stress by modulating the balance of phosphorylated/dephosphorylated species of the sphingolipid biosynthesis pathway.21 These enzymes are particularly intriguing due to their evident association with the aforementioned processes. The specific dynamics or interactions of de novo synthesis path that might involve these kinases, or how the enzymes' functions might impact broader cellular or physiological responses, has not been looked into. However, both free LCBs and phosphorylated forms typically occur in low abundance in plants, making it difficult to measure them.22

Considering the regulations and broad impacts of sphingolipid metabolism in diverse biological contexts, our study conducts a quantitative dynamic metabolic flux analysis on the de novo sphingolipid biosynthesis in Arabidopsis. While numerous investigations have delved into the regulation of the de novo sphingolipid biosynthesis, employing tools like isotope dynamic labeling in yeast and mammalian cells, the detailed dynamic metabolic flux analysis of de novo synthesis of sphingolipids in Arabidopsis remains relatively unexplored.1,23,24,25,26 As an example, Chen et al.27 developed a model rooted in pulse labeling, but it did not focus on de novo synthesis. This line of inquiry is essential, given the potential hazards of alterations in de novo sphingolipid biosynthesis. Well-documented cases, such as the effects of fumonisins—a mycotoxin found in maize—highlight these hazards, with outcomes ranging from diseases in plants to cancer in mammals.

To elucidate the sphingolipid pathway, stable isotope labeling has become an invaluable technique,24,28 leveraging tools such as mass spectrometry (MS) and nuclear magnetic resonance spectroscopy (NMR).29 Although 13C remains the popular choice for labeling, 15N and 18O isotopes are gaining momentum.29 In autotrophs, isotopically labeled nutrients, often introduced as nitrate salts, provide invaluable tracers for monitoring metabolic flows under real-life conditions.30,31,32 This method’s efficacy is highlighted in its application to organisms like Saccharomyces cerevisiae (baker’s yeast) and Arabidopsis, unearthing insights into lipid and photosynthetic processes. However, with sphingolipids, given their essential role in various cellular processes, accurate detection and quantification is necessary. Owing to their unparalleled sensitivities and specificities, MS tandem MS (MS/MS)-based approach has become a gold standard for detection and quantification tasks, ensuring precise measurements of these vital molecules in a variety of contexts. This analytical method’s precision not only allows for the clear identification of sphingolipids but also enables detailed insights into their intricate structures and molecular dynamics.28,30,31,33

In the multifaceted realm of sphingolipid metabolism and the dynamic nature of plant metabolism, tools for accurate exploration of these pathways are vital. Computational systems biology addresses this, either affirming established hypotheses or generating new ones from observed patterns. In recent times, constraint-based genome-scale metabolic models (cGMMs) have been crafted to probe the fundamental metabolism of several plant models, with Arabidopsis being a notable subject.34 Isotopic-labeling-based metabolic flux analysis (MFA) offers a quantitative, mechanistic portrayal of cellular phenotypes, producing a flux map grounded on the stoichiometric balance of isotopic labeling patterns.35 Typical steady state MFA offers an initial glimpse into labeled product distributions, but it does not directly discern intermediate turnover rates. This limitation is addressed by employing a detailed non-stationary isotope dynamic metabolic profiling, which mechanistically yields turnover rates and transient metabolic fluxes.29,36,37,38 A pioneering technique in this area is dynamic metabolic flux analysis (DMFA), which can be viewed as an enhancement of MFA,39 specifically, isotope tracer experiments incorporating time-dependent features for an exhaustive biological analysis. Within the realm of sphingolipid metabolism, utilizing DMFA to track temporal changes offers insights into metabolic stability, which is crucial for therapeutic explorations. A nuanced analysis of metabolic flux shifts can pinpoint specific molecular subspecies or enzymes, enabling targeted modifications to regulate downstream processes.

DMFA consists of two primary strategies: the static optimization-based approach (SOA) and the dynamic optimization-based approach (DOA) as outlined by Antoniewicz.39 SOAs simply average data, segmenting culture time into metabolic phases, or use data smoothing techniques for a stepwise solution approach. DOAs, on the other hand, use piecewise linear flux functions across continuous time segments. While comprehensive dynamic optimization would be ideal for accuracy, it is often impractical for extensive models. Conversely, SOAs may compromise on accuracy and solution stability. However, SOA has been employed in various biological organisms modeling, from Escherichia coli (E. coli) to Arabidopsis. DOA, with its distinct benefits, has also been used for E. coli metabolism and signaling pathways in S. cerevisiae.35,40 As exploration in this domain widens, both SOA and DOA have set the stage for groundbreaking studies on plant metabolism, such as resource allocation in Arabidopsis. Nonetheless, DMFA’s effectiveness is tethered to technological advancements. With the emergence of high-throughput metabolomics and computational modeling, dynamic models of metabolism (dMMs) have been established.41 Similar to DMFA, dMMs are also generally crucial for grasping the metabolic regulation of cellular perturbations over time durations. However, predicting the repercussions of such disruptions via dMMs remains crucial but challenging, considering the extensive regulatory and feedback mechanisms involved.

Moreover, with experimental metabolomics, particularly within the context of time-dependent tracer experiments, the inherent variability of data poses a significant challenge to accurately detecting molecules such as sphingolipids, bioactive lipids crucial for understanding cellular fate regulation. This variability stems from a range of sources, including biological variation, sample handling, and analytical processes.42 The task is further complicated in plants including Arabidopsis, which exhibit a higher degree of heterogeneity in sphingolipid LCBs compared to yeast and animals, making the identification of individual lipid species challenging due to identical mass spectrometry profiles and isotope interference.22,43 Although current sphingolipidomic studies primarily utilize liquid chromatography-tandem mass spectrometry (LC-MS/MS) alongside nuclear magnetic resonance (NMR) for compound identification and the creation of a multiple reaction monitoring (MRM) library, this approach faces limitations, especially in accurately measuring low-abundance sphingolipids due to significant background effects and variability.44 Despite these challenges, ongoing research efforts aim to develop an integrated approach that combines high-throughput analysis, precise identification, and extensive coverage, as evidenced by some studies reporting sphingolipid measurements across various cell types and organisms with coefficient of variation (CV) ranges from approximately 10%–40%.45,46,47

Addressing these challenges requires comprehensive analytical strategies and robust tools to ensure accurate and reliable data interpretation. The standard DMFA protocol is also prone to complications with significant variation in the data. To counteract this, methodologies incorporating regularization, ensemble modeling, and bootstrapping are being established to fortify the confidence in predictions derived from the DMFA approach. Such techniques are instrumental in handling the variations and uncertainties inherent in dynamic systems, paving the way for more reliable metabolic flux estimations. This advancement is exemplified by the work of Hebing et al.,48 which demonstrated the efficacy of DMFA algorithms that capitalized on regularization and constraint-based methods to address data variability, thereby obtaining more robust elementary mode (EM) flux predictions as indicated by wider 95% confidence intervals in their parameter fits. Our study further extends these advancements by focusing on the dynamics of intracellular metabolic pathways, employing targeted metabolomics and dynamic labeling techniques to provide a more nuanced understanding of metabolic fluxes at the intracellular level. This approach, focusing beyond the conventional analysis of uptake or secretion rates, offers a complementary perspective. Specifically, our methodology emphasizes the critical role of detailed intracellular analysis within the metabolic flux analysis toolkit. By directly tackling the challenges posed by temporal-intracellular noisy concentration measurements and the complexities of metabolic networks, it significantly advances our capacity to quantitatively analyze and understand metabolic systems.

In this work, we used a time-course 15N isotopic labeling technique coupled with targeted metabolomics to gain a dynamic and quantitative understanding of the de novo sphingolipid biosynthesis in T87 cell cultures. We delved deeper using a DMFA framework, which incorporated both regularization and dynamic flux sampling, termed r-DMFA. This approach allowed us to observe how enzyme activities change over time. Despite the challenge of dealing with varying data, our findings offer valuable insights into the behavior of sphingolipids within cells, an effort demonstrating the effectiveness of dynamic modeling in revealing the complexities of biological systems. Figure 1 provides a comprehensive view of our study workflow and methodology.Figure 1 Overview of the workflow and methodology used in this work

(A) Experiment design and labeling workflow depicting cell culture, growth, isotope labeling, quantification, and initial metabolomics analysis.

(B) Schematic illustration of the step-by-step procedure employed to obtain transient metabolic fluxes and analyze the inherent dynamic modulations of the Arabidopsis sphingolipid biosynthesis pathway.

(C) In silico knockout (KO)/essentiality hypothesis testing using the customized MATLAB implementation of the r-DMFA modeling framework.

Our computational analysis emphasized the pivotal role of SBH in preserving cellular equilibrium. A shortage of SBH led to distinct metabolic shifts and an increased predisposition to PCD, a consequence of sphinganine accumulation. Notably, active LCBK channel, which comprises its analogous and putative reversible enzyme, was identified to regulate shortage SBH activity and by potentially acting as a rate-determining factor in the degradation of phosphorylated phytosphinganine/phytosphingosine. Additionally, our in-depth flux analysis also showcased a cooperative relationship among FA2H, CSII, and GCS enzymes—key pillars of cellular stability. Disturbances in their functions jeopardized cell viability and raised PCD risks, as evident from potential LCB/Cer build-up and deviations from standard wild-type (WT) dynamics. From these insights, we suggest particular control nodes in the sphingolipid synthesis pathway. Our findings aim to create a predictive model of the sphingolipid network’s response to genetic and environmental changes, setting the stage for bioengineering stress-resilient crops.

Results

Comprehensive 15N profiling and isotopologue analysis of sphingolipid species

In the pursuit of a comprehensive method to quantify diverse sphingolipids in raw samples, mass spectrometry (MS) was chosen for its adaptability to individual analytes. Due to the complexity of plant sphingolipid structures, we combined MS with high-performance liquid chromatography (HPLC). This combination facilitated efficient transient sphingolipid measurements and conformed to the principles of the HPLC/electrospray ionization tandem mass spectrometry (HPLC/ESI-MS/MS) method as detailed in refs.22,43

In an in vitro experiment, unlabeled A. thaliana T87 cell suspension cultures established from the ecotype Columbia were grown in NT-1 liquid medium with the naturally occurring 14N isotope. 15N-labeled cultures were fed with NT-1 liquid medium containing K15NO3 and 15NH415NO3 (Cambridge Isotope Labs) to introduce a stable 15N isotope under the same conditions. HPLC/ESI-MS/MS was used to quantify 168 sphingolipid species over a time course. The incorporation of 15N into sphingolipids was monitored as a +1 m/z mass increase for both precursor and product ions for multiple reaction monitoring (MRM) transitions.

Following HPLC/ESI-MS/MS analysis, sphingolipid molecular species were measured until the mid-exponential growth phase, which occurs on day 6 (see also Figure S1). Figure 2B displays the total detected sphingolipids and highlights the abundance of the 14N and 15N isotopes for each species on day 5. This day is representative of the mid-exponential growth phase. The fingerprint patterns observed for each sphingolipid class in Figure 2B are consistent with the findings of Markham and Jaworski22 for crude samples of Arabidopsis thaliana Col-0, specifically of leaf tissue.Figure 2 Comprehensive 15N profiling and isotopologue analysis of sphingolipid species

(A) Overview of the sphingolipid biosynthesis pathway. The illustrated pathway model,49 tailored for r-DMFA analysis, captures the full spectrum of active metabolic transformations in Arabidopsis sphingolipid biosynthesis and acknowledges the dominance of either C16 or C24 fatty acids (FAs) in sphingolipids. Blue elements relate to C16 fatty acids, whereas red pertains to C24, with accompanying enzyme annotations (letters a–h) for specificity. SPT, serine palmitoyltransferase; KSR, 3-ketosphinganine reductase; SBH, LCB C-4 hydroxylase; CSI, class I ceramide synthase; CSII, class II ceramide synthase; SLD, LCB Δ8 desaturase; FA2H, fatty acid hydroxylase; GCS, glucosylceramide synthase; GlcCer, glucosylceramide; GIPCS, glycosylinositol phosphoceramide synthase; LCB, long-chain base; LCBK, long-chain-base kinase; DPL1, dihydrosphinganine phosphate lyase. Letters a–h indicate the enzyme reaction specificity to the leading reactive species. See also Tables S1 and S2.

(B) A quantitative snapshot of 15N-labeled sphingolipids in Arabidopsis T87 cell cultures. Saturated forms of dihydroxy and trihydroxy LCB groups are represented by “d18:0” and “t18:0”, whereas the unsaturated forms are represented as “d18:1” and “t18:1”, respectively. FA groups in saturated and unsaturated forms, for example with 16 carbon chain, C16, are represented as “c16:0” and “c16:1”, whereas the hydroxylated forms are “h16:0” and “h16:1”.

(C) Transient isotopologue fractions [(m + 1)/(m + 1 + m + 0)] of 36 sphingolipid species during 15N labeling (n = 4, two replicates per biological duplicates). Each depicted fraction represents the proportion of 15N enrichment in the measured analyte relative to its entire pool.

(D) Temporal abundance of 15N in sphingolipid species. Each profile showcases the corresponding temporal 15N abundances for the 36 sphingolipid species, to emphasize the de novo synthesis rate across biological replicates, A, B, C, and D. Cultures C and D are 15N-labeled replicates, whereas A and B are 14N unlabeled replicates. Data are represented as mean ± SEM.

To dynamically trace 15N incorporation into sphingolipid species, isotopologues, we analyzed the base m/z (m + 0) and +1 m/z shift (m + 1), of labeled and unlabeled sphingolipid species (hereafter termed “analyte”). The (m + 0) isotopologue, exemplified by 302.3 m/z for “d18:0”, represents the molecule’s mass where all atoms exist in their predominant isotopic form, notably 14N. In contrast, the (m + 1) isotopologue, such as 303.3 m/z for “d18:0”, denotes the mass shift resulting from 15N integration or from the natural presence of 13C, 18O, and 2H isotopes.

The relatively similar trends observed in the entire active pool measurements (m + 1 + m + 0) of most species across biological replicates under both labeled and unlabeled conditions (as depicted in Figure S2 and supported by Figure 2D) attest to effectiveness of this approach. Through dynamic 15N labeling, isotopologue fractions [(m + 1)/(m + 1 + m + 0)] were computed, which signify the proportion of 15N enrichment in the analyte relative to its entire pool. As such, isotopologue fractions during 15N labeling are internally adjusted to account for the quantity of pre-existing analytes prior to 15N integration, as depicted in Figure 2C. In more specific terms, the (m + 0) and (m + 1) isotopologues are equivalent measures for the abundance of N14 and N15 isotopes, respectively. Moreover, Figures 2C and 2D demonstrate the changing proportions of isotopologues for 36 different sphingolipid types, utilizing a 14N-labeled control to ensure the accuracy of our measurements across each sphingolipid species. Specifically, we established four biological replicates, A, B, C, and D. Cultures C and D are 15N-labeled replicates, whereas A and B are 14N-unlabeled replicates. We expected these controls to show consistent results, confirming the specificity of our measurements. The varying amounts of 15N observed over time, detailed in Figure 2D and Data S1, are indicative of de novo synthesis, which suggests active sphingolipid production from scratch. The LCB sphingolipids, critical for their limiting role in controlling the biosynthesis pathway’s pace upstream, were found in significant quantities.

The amounts (abundances, see Figure 2D) of phosphorylated LCBs (LCBPs) showed fluctuations, which could be due to their rare presence, inactive states, or the sensitivity of our measuring instruments. For sphingolipids marked by trihydroxy LCB structural backbone (“t18:0”/“t18:1”), we noted an increase in ions marked by 15N (again, reflected as m + 1) and a corresponding decrease in 14N (m + 0) ions, particularly noticeable from days 2–5 of the 15N labeling period. This pattern, similar to that seen in sphingolipids with dihydroxy LCB structures (“d18:0”/“d18:1”), confirms consistent labeling trends across different sphingolipid families.

Notably, there is also very high variability in measured levels of sphinganine (d18:0 LCB) and phytosphingosine (t18:0 LCB) for the majority of the time points (Figures 2D and S2). We further examined the temporal abundance profile for d18:0 LCB (Figure 2D, top left panel), which showed neither systematic nor instrumental noise but rather slight culture variation, particularly in culture C and D. This contrasted with the t18:0 LCB profile, for which we infer a random rather than systematic noise trend, especially for the labeled cultures C and D. This trend was supported by similar behavior in the total abundance of nitrogen (14N + 15N) in control cultures (A and B), suggesting transient fluctuations in the measurement instruments or environmental conditions that are not perfectly controlled. Examining the maximum coefficients of variation (CVs) in the intercultural assays (C and D) for the t18:0 LCB profile from day 2, excluding the lag phase, reveals significant biological variation between samples, with maximum CVs of 70% and 30%, respectively. However, the trends in variation, as presented in Figure S6, show that for the majority of species, CVs are relatively stable and stay below 20% across different time points. This suggests consistent measurements or the presence of systematic, rather than random, noise for these species. There are exceptions, such as LCBPs, t18:0 LCB, d180-hC24, d181-C16, GIP-d181-hC24, and GIP-d180-hC24, which display more variability and are present in lower abundance compared to other species.

Apparent turnover frequency distribution across sphingolipid species

In exploring the dynamics of de novo sphingolipid biosynthesis and turnover, we leveraged the metabolic pathway established in one of our prior studies49 (see Figure 2A; Table S1). Utilizing the in vitro 15N-labeled quantification derived from HPLC/ESI-MS/MS analysis, we assessed the metabolic turnover during the experiment by comparing the abundance of 15N-labeled sphingolipids to the total sphingolipid pool. The resulting dynamic isotopologue fraction (or isotopic incorporation/enrichment) curves (e.g., species in Figure 2C) demonstrate that labeled sphinganine (“d18:0”), associated with the restrictive uptake of serine, is integrated into sphingolipid before the exponential phase (specifically in the initial 24 h). The sphingolipid then displays relatively consistent turnover rates. The employed Arabidopsis T87 cell culture, which is undifferentiated and photoautotrophic, possesses a short doubling time of about 2.7 days. Notably, barring a few exceptions, sphingolipids with dihydroxy LCB backbones exhibit isotopic incorporation curves that gradually stabilize, reaching a plateau around day 5 of the exponential phase.

We propose an indirect method suggesting that the in vivo turnover number of a metabolic reaction can be apparently determined by observing the rate and extent of labeled nitrogen synthesis or integration. This approach provides a practical avenue for assessing catalytic speed and efficiency in living systems, particularly when direct measurements pose challenges. By applying first-order kinetic growth/decay models and regression analysis, as detailed in our STAR Methods section, we aim to ascertain the relative turnover numbers of de novo synthesized species. Hence, we calculated turnover frequency (TOF) distributions, also denoted as kcat+ (or simply k), by applying linear regression to the nonlinear model curve described in the STAR Methods section, [fraction(t)=1–e−kcat+t] yielding the distributions in Figure 3.Figure 3 Apparent turnover frequency (TOF) distribution across sphingolipid species

This figure illustrates the derived apparent turnover frequencies (TOFs or kcat+ values) for various sphingolipid species. Rapid turnover, relative to other species, is evident in GlcCer and GIPC species, indicating their pivotal role in dynamic cellular processes. Distributions are determined for species in biological duplicates (representative cultures C and D).

From the resulting distribution of catalytic efficiency rates kcat+, FA2Ha, responsible for converting d18:0-C16 to its corresponding hydroxy-ceramide (hCer) variant, demonstrated increased activity, as shown in Figures 2B and 3. A preference for converting the C16 fatty acyl to the C24 fatty acyl group was noted in LCB-Δ8-desaturase (see SLDb in Figure 3) activity. Class I ceramide synthases (CSIa&b) showed similar activities, suggesting a competitive dynamic among pathways, with CSII displaying a preference for C24 fatty acyl groups, particularly CSIIa. Moderate consistency in TOF distributions was observed across biological replicates, with a Bhattacharyya coefficient50,51 overlap average of 51.89% ± 30.10% across reactions. In the context of de novo synthesis, distinct clusters of turnover frequency (TOF or kcat+) values were identified among specific sphingolipid pairs, particularly between GCSf/h and GIPCSf/h. These pairs convert hCer t18:1-hC16/C24 into GlcCer t18:1-hC16/C24 and GIPC t18:1-hC16/C24, respectively, highlighting the metabolic demands on t18:1-hC16/C24 for maintaining cellular homeostasis. However, increased accumulation of hCer d18:1-hC24 substrate was observed, attributed to GCSg activity, acting as a potential regulatory shunt (Figure 2D). It is important to note that the cell line used lacks essential hormones for cellular differentiation, aligning with research on GCS’s role in cell-type differentiation in A. thaliana,13 thus underscoring the need for further research on the GCSg reaction.

Certain sphingolipid species, such as GlcCer and GIPC, undergo rapid turnover, indicating their involvement in dynamic cellular processes, whereas other species showed relative stability. Specifically, high turnover rates for GIPC species, especially those containing trihydroxy LCBs (t18:0 and t18:1), as shown by the elevated kcat+ values for GIPCSf and GIPCSh, highlight GIPC’s critical role in lipid metabolism and membrane stability in Arabidopsis.52 This observation is consistent with the known significance of tri-hydroxylated sphingolipids in plants, particularly the monounsaturated and saturated LCB derivatives (t18:1 and t18:0).3,52 Notably, the ratios of de novo pools between t18:1-hC24 and t18:0-hC24 for both hCer and GIPC maintained a consistent ∼1.5-fold difference, emphasizing the importance of these molecules in cellular balance. This finding also underscores the indispensable function of sphingolipid long-chain-base desaturases (SLD).

Transient de novo fluxes within the sphingolipid pathway network

Analysis of reaction fluxes over time revealed key processes crucial for growth, proliferation, and cellular homeostasis. Our study employed an enhanced r-DMFA modeling, incorporating cubic spline interpolation and Hessian matrix regularization, to effectively dissect dynamic shifts in metabolic fluxes and bolster the robustness and flexibility of r-DMFA fitting (see Figure 1 for workflow and see STAR Methods). Transient fluxes were determined using 15N-labeled quantification from HPLC/ESI-MS/MS analysis as a relative measure of de novo sphingolipid biosynthesis. The metabolic pathway model adopted for this DMFA analysis aligns with the one described in our previous work (Alsiyabi et al.49) (see Figure 2A; Tables S1 and S2). Similar to the referenced work, our model encompasses a thorough knowledgebase of all active metabolic transformations in Arabidopsis sphingolipid biosynthesis, accounting for enzyme promiscuity, parsimonious principle of minimal enzyme activity in most tissues, and the dominance of either C16 or C24 fatty acids in sphingolipids.11,53 We also adopted three specific assumptions from the referenced study, including the following: (1) LCB C4 hydroxylase (sbh gene) primarily hydroxylates free long-chain bases, (2) fatty acid α-hydroxylase (fa2h gene) targets only ceramide-bound fatty acids, and (3) LCB Δ8 desaturase (sld gene) produces unsaturated LCBs bound to ceramides. By adopting these assumptions based on previous research, we minimize the degrees of freedom. Additionally, akin to the referenced study, our model consolidates certain metabolic pathway steps due to a lack of relative kinetics information and excludes reactions converting non-hydroxylated ceramides into complex forms, consistent with the composition of most glycosylated sphingolipids.

Utilizing the r-DMFA framework, transient de novo fluxes of internal reactions in the network were predicted and are illustrated in Figure 4. These values were normalized to serine palmitoyl transferase (SPT) reaction’s flux to demonstrate how the flux through this reaction branches into different sphingolipid products. Metabolite pool sizes were also determined as an output of the DMFA fitting procedure (see STAR Methods). Interestingly, the majority of the produced long-chain bases (LCBs), approximately 93%, are directed through class II ceramide synthase to generate FA containing ceramides, particularly those with more of C24 (very long-chain fatty acid [VLCFA]) and t18:0-LCB.Figure 4 Transient de novo fluxes within the sphingolipid pathway network

(A) Metabolic network map used for r-DMFA analysis, highlighting internal reactions of the sphingolipid pathway upstream of the hydroxyceramide (hCer) metabolite junction.

(B) Corresponding transient fluxes from the pathway shown in (A).

(C) Transient fluxes for processes occurring after the hCer junction. Flux values were normalized to serine palmitoyltransferase (SPT) activity. A predominant channeling of long-chain bases (LCBs) (∼93%) through class II ceramide synthase for the production of very long-chain fatty acid (VLCFA)-containing ceramides, particularly with C24 and t18:0-LCB, is observed. Fluxes shown in (B) and (C) are the mean distributions for cultures C and D (n = 2). Data are represented as mean ± SEM.

The predictions indicate that the majority of the de novo species undergo recycling only after reaching the mid-exponential phase, as evidenced by the negative fluxes observed subsequently. The decision not to recycle these de novo species until the mid-exponential phase suggests a strategic metabolic shift that may be essential for cellular functions, potentially signaling a switch from biosynthesis to recycling that aligns with cellular demands for homeostasis and regulatory compliance. This could have profound implications for our understanding of sphingolipid dynamics in health and disease, where such regulatory mechanisms may be pivotal.

We note that our models induce a dynamic loop involving the LCB kinases, characterized by positive fluxes for LCBKa and negative fluxes (degradation) for LCBKb (refer to Figure 4B). Despite the presence of this loop, the models highlight the crucial function of SBH, which diverts approximately 80% of the metabolic flux away from the loop and largely unaffected by the minimal fluctuations observed in LCBP metabolomics data. Although we could have simplified the models by assuming LCB kinases (LCBKa/b) function in a unidirectional manner, we intentionally preserved their reversible activity. This choice was aimed at broadening the scope of our predictions, facilitating a more extensive investigation of metabolic pathways without being constrained to specific, expected outcomes. By adopting this flexible approach, we seek to enhance our comprehension of metabolic processes, steering clear of the restrictions posed by the presumption of irreversible enzyme actions. However, the variations noted earlier in the levels of LCB phosphates (LCBPs) are deemed relatively minor, especially when considered against a backdrop of changes less than one unit in magnitude compared to other metabolic profiles.

Analysis of enzyme demand and relative cost distribution in the sphingolipid biosynthesis pathway

To delve deeper into the enzyme requirements needed by the cell to match the experimental observation and maintain the observed WT flux distribution, the enzyme demand was calculated as a function of metabolite levels and transient de novo synthesis fluxes (v).54 This calculation provided a valuable measure of understanding resource allocation. The relative enzyme cost distribution computed is depicted in Figure 5A. This analysis clearly revealed a significantly low level of CS1 (LOH2), and GCS enzyme level, suggesting a close regulation of the enzyme’s function in response to substrate availability and metabolic demands, which could result in reduced specific activity. However, it is plausible that other enzymes in the pathway may compensate for this deficiency, ensuring the overall efficiency of the biochemical process.55 Furthermore, this relatively low enzyme activity warrants additional investigation to determine if environmental stressors or nutrient limitations could instigate higher demands, thereby linking these enzyme activities to an adaptive response aimed at conserving energy and resources.Figure 5 Analysis of enzyme demand and relative cost distribution in the sphingolipid biosynthesis pathway

(A) Distribution of relative enzyme costs, showcasing the energy and resource allocation within the metabolic network.

(B) Enzyme specificity investigation, representing the demand of each enzyme concerning each reactive sphingolipid species, highlighting the interplay of enzymes with various lipid species.

As depicted in Figure 5A, the enzyme cost for most of the biochemical transformations associated with biosynthesis (de novo) is relatively at par with the first and rate-limiting step of the pathway, the reaction catalyzed by SPT.4 This observation underscores the essential roles of SPT, SBH, and CSII (LOH1/3) in the thermodynamic feasibility of de novo synthesis. It also suggests that LCB branching at SBH may play a pivotal role in maintaining homeostasis within the sphingolipid pathway.55 Moreover, this bottleneck in enzyme demands, which can directly infer thermodynamic bottlenecks, indicates these enzymes are key players required to maintain the measured reaction rates and/or cellular homeostasis.

Nonetheless, an almost 2-fold relative cost was observed for FA2H, which is connected to the low turnover numbers observed in Figure 3. This wide range of distribution for FA2H also suggests that reactions involving this enzyme have a broad range of positive substrate/product concentration ratios, which can result in an overall forward flux of the reaction.55 This implies that a closely regulated mid-stream enzyme like FA2H may be a key player in de novo synthesis. Subsequently, a two-sample Kolmogorov-Smirnov test (KS test)56 was performed to determine if there were any statistically significant differences in the distributions of the calculated enzyme costs computed for experimental biological replicates. Interestingly, a significant difference (p < 0.001) was found for four enzyme classes, namely CS1 (LOH2), GCS, GIPCS, and SLD (t18:1). Visual inspection of these distributions (can be inferred from Figure 5) revealed that one replicate is a subset of the other, with calculated average Bhattacharyya coefficient overlap of 95.90% ± 3.91%.

Further examination of enzyme specificity was conducted by computing each enzyme demand concerning each reactive sphingolipid species, as shown in Figure 5B. This analysis revealed enzyme cost specificity with the highest visible variability in the FA2H class of enzymes. Notably high demands were observed with FA2Hd and FA2Hh, which convert t18:0-C24 and t18:1-C24 to their respective analogous hCers. This observation underscores the intrinsic role of FA2H in meeting the heightened metabolic demands of the t18:0/t18:1-hC24 species. Additionally, in general, enzyme specificities toward C24 VLCFA were observed to be higher than those toward FA C16. The average Bhattacharyya coefficient overlap for the biological replicates across the reactive species was 90.21% ± 7.80%. Interestingly, the low demand of GCS enzyme agrees with previous studies’ postulates linking GCS to cellular differentiation, which our cell-line lacks.13

Insights into metabolic pathway regulation through flux modulation ratio analysis

Building on these observations, a uniform convex polytope sampling strategy was employed to thoroughly analyze the high-dimensional dynamic flux solution spaces, considering different alternate optimal dynamic flux distributions. This involved uniformly drawing representative random transient concentration profiles within the given experimental constraints to truncate the bounded flux polytope, a method commonly used for characterizing the solution spaces of metabolism48,57,58,59 (see STAR Methods). The robustness of this approach is underscored by the strong Spearman’s rho (ρ) correlation coefficients exceeding 0.6, as detailed in Figure S4, which validates the predictive models against the sampled data across various species in C and D cultures. In addition, the residuals analysis showed a tendency to cluster around zero, indicating a good fit for the species within both C and D cultures (see Figure S3 for details).

Further, we examined a hypothesis, involving the partitioning of the sampled dynamic flux distribution set within the metabolic network. This hypothesis suggests that the relative flux modulation ratios, defined as the fold changes of each flux measurement relative to the mean flux of the distribution, directly influence the biological partitioning of reactions. Specifically, it was hypothesized that reactions with high variations in the distribution of flux modulation ratios would display a distinct partitioning pattern compared to those with low variations in the distribution of modulation ratios, thereby indicating a potential regulatory mechanism in the distribution of metabolic flux.

Moreso, cellular homeostasis, which is essential for the proper functioning of cells, is linked to the tight regulation of metabolic pathways. Figures 6A and 6B provide insights into how variations in flux modulation ratios across different biological replicates could be linked to the maintenance of cellular homeostasis. Specifically, the red cluster/partition, characterized by low variations in flux modulation ratios, encompasses reactions that are tightly regulated and likely constitute the limiting set of the biosynthesis (de novo). The identified candidate reactions (CSIIb, SLDd, FA2Hh, and GIPCSh) are in line with the activities of known essential enzymes such as SPT, KSR, and SBH. This tight regulation is crucial for maintaining cellular homeostasis, as it ensures that the synthesis of key cellular components is carried out efficiently and consistently, supposedly under varying environmental conditions. On the other hand, the green and blue clusters/partitions, characterized by moderate to high variations in flux modulation ratios, consist of metabolic invariant sets that have little to no tight regulations on de novo synthesis (see Figures 6C and 6D boxplot distributions). These reactions, while not as tightly regulated, still play a vital role in cellular homeostasis by providing flexibility and adaptability to the metabolic network, allowing cells to respond to changes in the environment while maintaining their internal equilibrium. Together, these findings highlight the critical balance between tight regulation and flexibility in the metabolic network, which is essential for maintaining cellular homeostasis and ensuring the proper functioning of cells.Figure 6 Insights into metabolic pathway regulation through flux modulation ratio analysis

(A) Scatterplot of flux modulation ratios for the first—culture C—biological replicate, depicting variations and highlighting tightly regulated (red cluster) and flexible metabolic reactions (green and blue clusters). Number of samples, n = 4E6.

(B) Scatterplot of flux modulation ratios for the second—culture D—biological replicate.

(C) Boxplot distribution of flux modulation ratios for the first—culture C—biological replicate.

(D) Boxplot distribution of flux modulation ratios for the second—culture D—biological replicate, providing a view of metabolic congruence or divergence.

It is important to note that the robustness of our predictive insights is not compromised by the fluctuations in the metabolomics data, as evidenced by the consistent clustering observed in the DBSCAN analysis across a broad range of sample sizes, detailed in Figure S5. This steadfastness amid varying sample volumes—from 1E+4 to 4E+6—reinforces the resilience and precision of our models in reflecting the true dynamics of metabolic pathways, further validating the predictive power and reliability of this approach.

Additionally, a customized MATLAB implementation of the r-DMFA modeling framework was developed to precisely predict transient metabolic reaction fluxes following enzymatic perturbations as knockout (KO) (see STAR Methods). This in silico KO/essentiality hypothesis testing approach could potentially minimize extensive experimental iterations in subsequent investigations and validations. This process was repeated for each enzyme class and reaction specificity and can be extended to probable enzymatic perturbations. As briefly depicted in Figure 1B, the experimentally fitted and transient metabolomics result in the predicted transient reaction fluxes for the reference (wild type). The obtained WT profiles are then used to generate dynamic models whose metabolic responses/profiles are screened with respect to different perturbations.

In silico metabolic analysis of single enzyme deletions and reaction knockout perturbations

To evaluate the metabolic network’s functionality without any additional regulatory interactions, we induced perturbations via single enzyme and reaction knockouts. We then leveraged the constraint-based r-DMFA framework, in similitude of established global minimization of metabolic adjustment (MOMA) principle,60 to predict each dynamic model’s response to the induced perturbation. As an objective here, we compared the predicted responses from single knockouts with experimental WT predictions, such that the r-DMFA framework, utilizing quadratic programming and relaxing optimal WT flux assumptions, more accurately predicts metabolic phenotypes of mutants and assists in comprehending cellular adaptation to enzymatic dysfunction. Our methodology, which utilizes the convex transformation methods inherent in the global optimization algorithm of the BARON solver,61 generates predictions that serve as a preliminary guide before conducting laboratory experiments.62

The metabolic network model was first examined through the lens of single enzyme deletions or knockout (KO) perturbations. By incorporating specific constraints into the model and applying an evaluation metric, zd, normalized to the maximum distance on a logarithmic scale (as depicted in Figure 7A and detailed in the STAR Methods section), the metabolic distance between each in silico mutant dynamic flux phenotype and the WT dynamic flux phenotype was assessed using Spearman’s rank correlation, as shown in Figure 7B. This approach was essential to understand the role of each enzyme within the metabolic network and to pinpoint those critical for the WT phenotype. The analysis in Figure 7A highlighted the significant roles of enzymes such as SPT, KSR, LOH1/3 (or CS II), FA2H, SLD, and GIPCS, which displayed pronounced metabolic distances from the WT, with zd values greater than 0.5-fold of the limiting SPT enzyme activity. LOH2 (or CS I) and GCS had moderate zd values around 0.5-fold.Figure 7 In silico metabolic analysis of single enzyme deletions and reaction knockouts

(A) Evaluation of the metabolic network model through the metric zd, normalized on a logarithmic scale, to understand the role of single enzyme deletions or knockout (KO) perturbations.

(B) Spearman’s rank correlation analysis assessing the metabolic distance between each in silico mutant and WT dynamic flux phenotypes, visualizing ratios of positive and negative correlations between mutant and WT profiles.

(C) Analysis of single reaction knockouts, showing metabolic distance analysis.

(D) Correlation analysis at the reaction level for single reaction knockouts.

Interestingly, the SBH mutant exhibited a minimal (<0.5-fold) metabolic distance to the WT, contradicting its traditionally understood role in cellular stability. Its deficiency is projected to lead to metabolic distant phenotypes and an increased potential for PCD due to sphinganine (d18:0) accumulation.63 This observation necessitated further examination of flux activities within these mutants. Previous research has shown that SBH mutants, which lack trihydroxy LCBs, exhibit stunted growth, an inability to switch from vegetative to reproductive growth, heightened expression of PCD genes, increased sphingolipid levels, and altered fatty acid composition.63

Further analysis revealed two potential pathways for compensating for the lack of trihydroxy LCBs: (1) a bypass involving LCBK, with positive LCBKa activity and a hypothetical reversible LCBKb activity to produce trihydroxy LCBs, and (2) a more complex pathway potentially involving DPL1 in a reversible reaction to synthesize trihydroxy LCB-1P. This pathway might include steps such as reductive amination, hydroxylation, and phosphorylation. The synthesis of the initial LCB from C16-fatty aldehyde and phosphoethanolamine is presumed to be relatively straightforward. However, the subsequent introduction of hydroxy groups and the final phosphorylation step are expected to be more complex, requiring precise control over the hydroxylation’s regioselectivity and stereoselectivity to successfully produce trihydroxy LCB-1P. These identified flux modes suggest a shift toward adaptive regulation rather than inherent regulation that does not depend on phosphorylated trihydroxy species. Assuming the prior synthesis of trihydroxy species, the putative sphingoid-phosphate phosphatase (LCBPP1) is hypothesized to dephosphorylate these molecules, converting them back to their LCB precursors.

Additionally, ratios of reactions showcasing both positive and negative correlations between the in silico mutant and WT profiles were calculated, as illustrated in Figure 7B. These data provided insights into metabolic congruence or divergence. The observations confirmed the prominent roles of the aforementioned limiting mutants, with significantly low ratios (<0.5). However, a stimulated metabolic divergence was observed for CS II or LOH1/3, GIPCS, and SLD mutants, evident in the relative increase of negatively correlated reactions when compared to the WT profiles. Moreso, the low ratios observed in Figure 7B for FA2H, GCS, CS I, and CS II as well as SLD mutants suggest significant metabolic reprogramming indicative of synergistic role these enzymes may play along this pathway.

Figure 7B also suggests the buffering activities of LCBK, DPL1, and SBH mutants having high ratio of significantly positive correlated reactions at par with the observations from their zd values in Figure 7A. Given the observed low ratio of negatively correlated reactions in the DPL1 mutant, it suggests that DPL1 might be involved in some metabolic reprogramming. This further points to a reversible action of the LCBK, potentially a putative sphingoid-phosphate phosphatase (LCBPP1), highlighting it feasible as a key modulator of SBH function. The LCBK catalyzes the phosphorylation of LCBs to form bioactive LCB phosphates. These phosphorylated species are also known to modulate cellular growth and differentiation.

The analysis of enzyme specificity regulations was extended to single reaction knockouts, as depicted in Figures 7C and 7D, which highlight both the metabolic distance and correlation analysis at the reaction level. This layer of analysis aimed to pinpoint critical reactions within the metabolic network that, when perturbed, induce significant deviations from the WT dynamic flux phenotype, akin to the analysis at the enzyme level. Trends analogous to those observed in the enzyme level knockout were uncovered. A notable observation was the reconfirmation of enzyme specificities; the specificities toward C24 VLCFA were found to be higher than those toward FA carbon length of 16, C16.

Moreover, the activity of SLDc bolstered the experimental observation that the concentration of VLCFA-containing ceramide d18:1-hC24 increases, thereby suggesting a reduction in the overall amount of LCBs produced due to the diminished rate through SPT. This reiteration underscores the interconnectedness of the enzyme abundance and specificity with respect to reactive species in a metabolic network. Understanding the regulations behind this balance will be key to engineer crops with improved yield while enabling robust response to biotic and abiotic stress. Reinforcing our observations on enzyme abundance and specificity, we rehypothesize that enzyme specificity plays a pivotal role in the rapid evolution and regulation of de novo synthesis, as evidenced by dynamic flux distribution in the metabolic network, guiding crucial metabolic regulatory shifts.

Discussion

Sphingolipids occupy a central position in both the biosynthesis and catabolism of cells and are highly dynamic during rapid growth.64,65 This dynamic nature projects large variations on the profiles of the species, especially the intermediates, including ceramides and LCBs, which are also bioactive molecules. Such variability is commonly observed in lipidomic studies, particularly those using MS-based measurement of low-abundant lipids. We have consistently observed this variation in our Arabidopsis cell line across several unpublished experiments.66 In this study, we develop and implement a regularized dynamic metabolic flux analysis (r-DMFA) framework in a bid to understand the dynamic nature of sphingolipid synthesis despite the challenges posed by experimental variability.

Isotopic labeling for sphingolipid profiling

We demonstrate the efficacy of 15N stable isotope labeling combined with HPLC/ESI-MS/MS for dynamic profiling of sphingolipid species in Arabidopsis thaliana. The use of 15N labeling allowed precise tracing of nitrogen incorporation into sphingolipids, revealing significant insights into sphingolipid metabolism during the mid-exponential growth phase. This method facilitated precise attribution of m+1 shifts to 15N incorporation, excluding contributions from natural isotopes. This approach offers a valuable tool for studying sphingolipid biosynthesis and turnover in plant cells. Each molecule of LCB, Cer, hCer, and GlcCer contains one nitrogen-15 atom. While Arabidopsis GIPCs may occasionally feature a second nitrogen atom on one of their sugar residues, this modification is uncommon in Arabidopsis leaves, and its frequency in cultured cells is unknown. The presence of a second 15N in certain GIPC varieties does not impact these results, but it should be noted.

A primary challenge in quantifying sphingolipids using 13C-labeling via MS involves issues in transitions, modifications, and fragmentation.1 To address this, 15N stable isotope labeling was utilized, targeting the nitrogen atom in the LCB, to provide a direct quantification approach for various sphingolipid species. The contribution of naturally occurring 13C isotopes is about 1%—meaning only 1 out of 100 naturally occurring carbon atoms is 13C—approximately 30% of molecules containing around 30 carbon atoms will include at least one 13C atom. For molecules with 18–48 carbon atoms, the probability of having one or more 13C atoms ranges from approximately 16.9% to 38.8%. This range suggests that for larger sphingolipids, the likelihood of 13C impact is significant. This is particularly relevant for larger sphingolipids such as GlcCers and GIPCs. This is evident in the relatively higher 15N (m+1) measurements for GIPCs in the control cultures (Figure 2D). However, having a 14N control measurement allows us to distinguish between the effects of naturally occurring isotopes and the specific incorporation of 15N. By comparing the isotopologue distributions in the 14N-labeled controls with those in the 15N-labeled samples, we can accurately attribute observed m+1 shifts to 15N incorporation rather than 13C or other natural isotopes. These control data are critical for interpreting the dynamics of 15N incorporation into sphingolipids accurately.

To address notable fluctuations within technical replicates for each biological sample, we implemented a tailored analytical strategy to accurately capture the significant intrasample variability present in our data, treating each biological replicate independently. This approach was particularly crucial due to the substantial variation highlighted by the incorporation ratio bands in Figure 2C, prompting us to utilize bootstrapping for our analyses. This statistical method allowed us to quantify the variability within the standard deviation of replicate measurements, effectively illuminating the inherent variability in our dataset. However, the reliability of the incorporation data for this substrate remains limited, advising caution in drawing broad conclusions about its regulatory impact.

Dynamic sphingolipid metabolism during rapid growth

Our findings highlight the critical role of sphingolipids in plant cellular processes such as growth, differentiation, and stress response.4,10 The mid-exponential growth phase, characterized by peak cell division and growth, showed notable changes in sphingolipid metabolism, reflecting increased demand for membrane components. Focusing on this phase provides a foundational dataset for further research on sphingolipid dynamics and their physiological implications. Despite observed variability, the approach’s robustness is evident in the consistent trends across biological replicates. The steady-state integration of 15N into sphingolipids by day 5, reaching approximately 80%–100% (Figures 2C and 3), underscores the method’s sensitivity and accuracy. In addition, the isotopic abundance of sphingolipids at mid-exponential growth on day 5 (Figure 2B) aligns with previous studies,22 validating the methodology and offering new insights into sphingolipid biosynthesis in A. thaliana.

Moreover, despite the inherent complexities of metabolic networks, our observations suggest a minimal influence from factors such as enzyme compartmentalization and substrate availability. This inference is supported by our analysis of 15N abundance relative to the entire nitrogen pool, which indicated negligible issues regarding substrate availability, with the notable exception of d18:1-hC24. The pattern of 15N accumulation and its incorporation across nearly all sphingolipid species suggests a widespread occurrence of de novo synthesis. Our modeling, based on the assumption of a first-order exponential relationship between the labeled nitrogen incorporation ratio and turnover rates, aligns with observed trends. The isotopic incorporation data (Figures 2C and 3) show an exponential increase in (m + 1) isotopologue fractions for most sphingolipids up to the mid-exponential phase (Figure 2C). Interestingly, not all sphingolipid active pools, specifically t18:1-hC24, follow an expected exponential growth pattern (Figure 2D). Instead, analysis shows a direct replacement of 14N with 15N, keeping the overall pool size trend relatively stable with unique variation (Figure S2). This suggests a dynamic balance of this essential species in Arabidopsis. Further research is needed to elucidate the dynamics and implications of this substitution effect within sphingolipid pools, particularly in understanding the metabolic turnover and stability of sphingolipids like t18:1-hC24.

Enzyme roles in the metabolic network regulation

The activity of FA2H in converting d18:0-C16 to hCer variants suggests a critical role in metabolic regulation, potentially involved in mediating PCD.6,67 This pattern, observed in a separate experiment with the same cell line, validates the observation and rules out contamination (unpublished data). In addition, a study by Ukawa et al.68 has emphasized the critical role of the hydroxylation step of sphingolipid fatty acids, which introduces essential hydroxyl groups for forming nanodomains within the plasma membrane. These nanodomains, primarily composed of sphingolipids and sterols, are fundamental to the plasma membrane’s structure, underscoring the significance of FA2H’s role in cellular organization.

Despite notable variability captured by the Bhattacharyya coefficient, the turnover frequency (kcat+), for LCB C4-hydrolase (SBH) and LCB Δ8 desaturase (SLDb)—presumed to synthesize the unsaturated LCBs backbone in the trihydroxy phytoceramides ceramides (t18:1-C16/C24) and not the free LCBs69—exhibited substantial overlaps of over 94% between replicates. This emphasizes SBH’s critical function in diversifying di- and trihydroxy LCBs. Our findings align with the known significance of tri-hydroxylated sphingolipids in plants and underscore the importance of sphingolipid LCB desaturases (SLD) in maintaining cellular equilibrium. In addition, the kcat+ value distribution for the CSII and GCS enzymes revealed significant differences between replicates, suggesting their differential and regulatory roles in the sphingolipid metabolism. However, the rapid turnover of specific sphingolipid species, such as GlcCer and GIPC, underscores their pivotal role in dynamic cellular processes, including signaling and membrane dynamics. These findings suggest that these sphingolipids are integral to the cell’s ability to adapt swiftly to environmental changes.

Furthermore, our in silico knockout simulations identified the differential activity of LCBK, which could mediate loss of SBH function by either serving as a recycle shunt or playing a limiting role in degrading phosphorylated phytosphinganine (LCB-1-phosphates [LCB-1Ps]). Previous studies5,18 have inconclusively explored the interaction between LCBK and its reversible counterpart. We hypothesize that the putative sphingoid-phosphate phosphatase (LCBPP1/SPP1), which dephosphorylates these molecules back to their LCB precursors, acts as a potential rate-limiting enzyme in the degradation of LCB-1Ps via the competing path of dihydrosphinganine-1-phosphate lyase (DPL1),70,71,72 suggesting coordinated regulation. More comprehensive research is needed to understand the interactions between these enzymes and their role in maintaining cellular homeostasis. LCB-1Ps are recognized as lipid messenger molecules implicated in plant stress signaling.16 A recent investigation by Nakagawa et al.73 explored the role of LCBP phosphatase in regulating LCBP levels and its participation in the stomatal response to abscisic acid (ABA) signaling pathways in Arabidopsis, with potential implications for enhancing water use efficiency and drought resistance. LCBKs also promote immunity against bacterial pathogens,74 further highlighting the dynamic balance of LCBs/LCB-1Ps.

Moreover, our model predicts enzymes like GCS, LOH (2, 1/3), SLD, and FA2H have varying impacts on the metabolic network, from low through moderate to significant. Notably, the in silico mutant for the GCS enzyme, despite its relatively low enzymatic cost in the WT phenotype, exerts a considerable reprogramming effect on the metabolic network. This indicates that even less costly enzymes in the WT can profoundly influence the metabolic network, highlighting the system’s complexity. The relationship between enzyme abundance and metabolic fluxes is intricate, influenced by factors beyond expression levels. Integrating data on enzymatic activities, regulatory mechanisms, and metabolite levels is essential to accurately model and understand metabolic fluxes. Hence, the dynamic metabolic trajectory of enzymes and reactions within the sphingolipid biosynthetic pathway in Arabidopsis offers a captivating look into cellular metabolism. Through meticulous in silico analyses, we can elucidate the regulatory mechanisms underpinning this pathway.

Overall, the metabolic turnover rates established the pivotal roles of several enzymes, including SPT, KSR, LOH1/3, and FA2H. Dynamic metabolic flux analysis (r-DMFA) showcased the branching nature of metabolic flux, emphasizing the significant roles of particular enzymes in the WT phenotype. The GCS enzyme, identified as tightly regulated with low enzyme costs, exerts a considerable reprogramming effect on the metabolic network, highlighting system complexity and intricate connections. Dynamic flux estimates suggested that LCBK-mediated SBH activity shortage might be a potential regulatory scheme. Additionally, SBH, LCBK, FA2H, CSII, and GCS play crucial roles in ensuring cellular viability and managing programmed cell death. Enzyme specificity toward phytosphingosine-based (tri-hydroxylated) groups and VLCFA emerged as critical, influencing the overall metabolic trajectory. High-throughput metabolomics and advanced in silico modeling present a promising avenue for generating experimentally testable hypotheses, illuminating sphingolipid biosynthesis and setting the stage for metabolic engineering to improve crop yield and stress response mechanisms.

Advancements, applications, and future directions

With advancements in high-throughput metabolomics and in silico modeling, dynamic models of metabolism can enhance our understanding of metabolic regulation beyond the sphingolipid biosynthesis pathway. The r-DMFA framework can be applied to other systems to study the combined regulatory impacts of biosynthesis and catabolism of metabolic reactive species over time. This dynamic modeling framework will produce more experimentally testable hypotheses regarding regulation in the sphingolipid biosynthesis pathway, thereby minimizing exhaustive experimental efforts.

In conclusion, we underscore the robustness of our methodology, specifically its adept handling of inter-sample variability—a formidable obstacle in metabolic studies that can mask the delicate dynamics of metabolic fluxes. By employing regularization, normalization, and ensemble techniques, our approach transcends the challenge of temporal variation within samples, enabling the extraction of significant patterns from complex datasets. Adjusting for variability both within and across samples has sharpened our ability to unveil the underlying biological mechanisms with greater precision and assurance. Our approach is substantiated by the marked improvement in the calculated average Bhattacharyya coefficient overlap of the distributions, which surged from 51.89% ± 30.10% to an impressive 95.90% ± 3.91% across the reactions, further illuminating the inner workings of sphingolipid biosynthesis and paving the way for future metabolic engineering endeavors.

Additionally, future research could explore the application of this approach to other plant species and developmental stages, as well as the potential impact of environmental factors on sphingolipid metabolism. Understanding the dynamic nature of sphingolipid synthesis will enhance our knowledge of plant physiology and stress responses, contributing to advances in agricultural biotechnology and crop improvement.

Limitations of the study

This study faces a few limitations that warrant attention. Firstly, it was conducted using cell suspension cultures, which may not fully replicate the conditions or behaviors of intact plants. Despite the inherent variability in experimental data, we have developed a protocol aimed at identifying common metabolic constraints amid noisy data, facilitating the study of their effects on metabolism. The challenge of obtaining noise-free, time-resolved measurements is particularly pronounced in the context of sphingolipid metabolism, as demonstrated by our initial turnover frequencies (TOF)/catalytic efficiency (kcat+) analysis, which is independent of the r-DMFA fitting technique used to explore metabolic solution spaces. We anticipate that sphingolipid synthesis exhibits significant sensitivity and variability in response to growth and cell signaling, a core focus of our investigation. By being able to first capture these early dynamics (prior to mid-exponential growth phase), we aim to deepen our understanding of how sphingolipid metabolism adapts across different growth phases and signaling contexts.

Furthermore, our computational (in silico) predictions of enzyme perturbations await experimental validation, which is unfortunately not within the scope of this study. We acknowledge that further experiments are essential to validate our predictions. Ideally, the least of these experiments would involve in vitro labeling studies using 2D-NMR for gene knockouts and generating temporal experimental metabolomics and proteomics data. However, these experiments are resource-intensive, time-consuming, and currently beyond the scope of this work. In future, such validation might also rely on cell suspension cultures, with subsequent studies potentially extending to intact plants through direct laboratory experiments. The extrapolation of data from cell suspension cultures to real plants poses a significant question for future research. Exploring the applicability of these findings to intact plants could illuminate pathways to engineering crops with enhanced stress resilience. Further investigation into this area is crucial for translating in vitro insights into tangible agricultural advancements.

STAR★Methods

Key resources table

REAGENT or RESOURCE	SOURCE	IDENTIFIER	
Chemicals, peptides, and recombinant proteins	
	
K15NO3	Cambridge Isotope Laboratories	Cat#NLM-765-1	
15NH415NO3	Cambridge Isotope Laboratories	Cat#NLM-390-1	
Murashige and Skoog	Phyto Tech Labs	M519	
Myo-inositol	Sigma-Aldrich	Cat#I7508	
2,4-Dichlorophenoxyacetic acid (2,4-D)	Research Products International	Cat#D43050	
Sphingolipid standards	Avanti Lipids -  Markham and Jaworski22		
	
Deposited data	
	
Data used for modeling	Zenodo repository75	Zenodo: https://doi.org/10.5281/zenodo.12811205	
	
Experimental models: Cell lines	
	
Arabidopsis thaliana T87 cells (Col-0)	Arabidopsis Biological Resource Center	Cat#CCL84839	
	
Software and algorithms	
	
Shimadzu Prominence UPLC system	Shimadzu Corporation	Model: LC-30AD	
4000 QTRAP mass spectrometer	AB SCIEX	Model: 4000 QTRAP	
Analyst 1.5 software	AB SCIEX	–	
MultiQuant 2.1 software	AB SCIEX	–	
MATLAB (r2021b)	MathWorks	https://www.mathworks.com/products/new_products/release2021b.html	
Custom kinetic model script	Developed in-house	Available at: https://github.com/ssbio/r-DMFA	
GitHub repository for r-DMFA implementation	GitHub	Available at: https://github.com/ssbio/r-DMFA	
Python version 3	Python Software Foundation	https://www.python.org	
CPLEX solver	IBM	https://ampl.com/products/solvers/solvers-we-sell/cplex/	
BARON Matlab-Interface	MINLP	https://minlp.com/matlab-baron-interface	

Resource availability

Lead contact

For additional information or to request resources, please contact Dr. Rajib Saha at rsaha2@unl.edu.

Materials availability

This study did not generate new unique reagents.

Data and code availability

• The data used for the r-DMFA study is publicly available on Zenodo at https://doi.org/10.5281/zenodo.12811205 and can also be found in the GitHub repository at https://github.com/ssbio/r-DMFA/tree/main/data.

• All original code, including those for statistical analyses, is available in the GitHub repository at https://github.com/ssbio/r-DMFA and can also be accessed via Zenodo at https://doi.org/10.5281/zenodo.12811205.

• The lead contact will provide any additional information needed to reanalyze the data reported in the paper.

Experimental model and subject details

Plant material and culture conditions

This study utilized undifferentiated and photoautotrophic Arabidopsis thaliana T87 cells originated from the Columbia ecotype, Col-0. Cells were cultured in NT-1 liquid medium (Murashige and Skoog medium enriched with vitamins, 30 g/L sucrose, 1 mM KH2PO4, 1 mg/L thiamine, 100 mg/L myo-inositol, and 2 μM 2,4-dichlorophenoxyacetic acid; pH 5.8 adjusted with KOH) in 100 mL flasks. Cultures were maintained at 22°C, 55% humidity, agitated at 120 rpm, and continuously illuminated at a photosynthetic photon flux density (PPFD) of 100 μmol m−2 s−1.

For internal standard and dynamic labeling, T87 cells were cultured in 15N-enriched NT-1 liquid medium containing K15NO3 and 15NH415NO3 (Cambridge Isotope Labs), under identical culture conditions. Standard sub-culturing was performed every 7 days, with 1 mL of the existing cell suspension transferred to 50 mL of fresh NT-1 medium. The 15N medium was introduced dynamically with daily sample collections initiating immediately and continuing for 6 days, up to the experimentally determined mid-exponential growth phase. A control culture was concurrently maintained in NT-1 liquid medium with the natural 14N isotope, under identical conditions. Each culture setup entailed two technical replicates with two biological replicates per setup.

To ascertain growth phases, growth rates were ascertained by freeze-drying and dry-weighing equal volumes of suspension cultures (bi-replicates) over a 15-day period. Data were then plotted on a semi-log graph, with linear regression applied to determine the growth rate from the slope of the best-fit line (see Figure S1 for growth curve).

Sphingolipids extraction and quantification

The method employed for sphingolipid extraction and quantification was adapted from Markham and Jaworski (2007).22 The profiling methods which have been modified was published in Cahoon et al. (2021).43 Specifically, 10–30 mg of lyophilized cells were homogenized and extracted employing a solvent mixture of isopropanol:heptane:water (55:20:25, v/v/v). Following extraction, the supernatants were dried and de-esterified using methylamine in a ethanol:water solution (70:30, v/v). The lipid extracts were then re-suspended in a solvent composition of tetrahydrofuran:methanol:water (5:2:5, v/v/v) with 0.1% (v/v) formic acid. Sphingolipid species analysis was conducted using a Shimadzu Prominence ultra-performance liquid chromatography system and a 4000 QTRAP mass spectrometer (AB SCIEX). Internal standards were incorporated for precise sphingolipid class identification.

For precise quantification of 14N/15N labeled sphingolipids, standard curves were generated through sphingolipid (analyte) isolation and analysis utilizing the 4000 QTRAP LC-MS/MS system across multiple dilutions. Specifically, the HPLC/electrospray ionization tandem mass spectrometry (HPLC/ESI-MS/MS) method, as described by Markham and Jaworski,22 was employed, predicated on the principles of ion transitioning, separation, and fragmentation. In this method, chromatographic separation was achieved using a reversed-phase C18 column, facilitated by a specified gradient elution of aqueous and organic mobile phases. Following separation, the eluent was subjected to electrospray ionization (ESI) to generate charged molecular ions, which were subsequently directed into the tandem mass spectrometer (MS/MS). Within the MS/MS, ions were initially filtered by the first quadrupole (Q1) based on their mass-to-charge ratios (m/z), then fragmented in the collision cell (Q2) to yield a spectrum of product ions, which were filtered and detected by the second quadrupole (Q3). The methodology hinged on the selection of optimal precursor-product ion transitions, accompanied by the meticulous optimization of mass spectrometer parameters including de-clustering potential (DP), collision energy (CE), and retention time.

Utilizing the HPLC/ESI-MS/MS technique, 168 distinct sphingolipid species/analytes (LCBs, ceramides, GlcCers, and GIPCs) were quantified over a designated time course. To dynamically trace 15N incorporation into sphingolipid species, analytes were scrutinized for base m/z (m + 0) and +1 m/z shift (m + 1) isotopologues. The (m + 0) isotopologue, e.g., 302.3 m/z for “d18:0”, denotes the mass wherein all atoms exist in their predominant isotopic form, predominantly 14N. In contrast, the (m + 1) isotopologue, e.g., 303.3 m/z for “d18:0”, represents the mass shift due to a single 15N integration. While the contribution of naturally occurring 13C isotopes is about 1% --- meaning only 1 out of 100 naturally occurring carbon molecules is 13C – approximately 30% of molecules containing around 30 carbon atoms will include at least one 13C atom. For molecules with 18–48 carbon atoms, the probability of having one or more 13C ranges from approximately 16-39%. This range suggests that the likelihood of 13C impact is significant, particularly for larger sphingolipids such as GIPCs, where the 13C isotope effect creates a notable proportion of (m+1) shifts. Despite this, the incorporation of 15N is pronounced and easily distinguished from possible 13C isotope effects using 14N-labeling control measurement, especially in the smaller sphingolipid classes that contain only 18 to 48 carbons. Additionally, technical duplicates were executed for each biological replicate inclusive of cultures grown on 15N and control cultures on 14N. Data analysis and quantification were performed using Analyst 1.5 and MultiQuant 2.1 software, with considerations for exact mass and retention time as described by Markham and Jaworski.22

Method details

Dynamic modeling

Dynamic metabolic modeling captures both the metabolic and regulatory state of a pathway through temporal reaction flux changes. Unlike many genome-scale metabolic models (GSMMs) which are static, dynamic metabolic models (DMMs) capture the dynamic nature of metabolic networks which is crucial for understanding and engineering biological systems. In this section, we highlight the in silico dynamic modeling approaches used to analyze the metabolic data and generate insights.

Estimating apparent turnover frequency (TOF) distributions

To estimate the Apparent TOF distributions, we employed a kinetic modeling approach that leverages the time-dependent incorporation of 15N-labeled nitrogen into sphingolipid species. The core of this approach is the application of a first-order kinetic model, which makes the following assumptions.(1) The rate of isotopic 15N incorporation into the species (rate of appearance and hence, turnover) equals the depletion rate of 14N from the pools of nitrogen (N) per species. Given the depletion rate of 14N is directly proportional to the 14N concentration at any given time. This follows the first-order irreversible kinetic relationship as: N14→kN15

∴d[N15]dt=−d[N14]dt=k∗[N14]=k∗([N]total−[N15])

Where, [N]total is the total concentration, N14+N15(2) The first-order rate constant k serves as a direct measure of catalytic efficiency, as it quantifies the rate at which a substrate is converted to product per enzyme molecule under saturating substrate conditions. Thus, the catalytic efficiency of the reaction enzymes is equivalent to the kinetic rate of isotopic labeling, reflected by the turnover number (kcat+≡k) in the reaction’s first-order rate law.

(3) At the start of labeling (t=0), the depleting species ratio, i.e., the ratio of unlabeled to total sphingolipids (N14N14+N15) is approximately one, assuming no N15 presence, while the ratio of labeled to total sphingolipids is zero. i.e., N14N14+N15(t=0)=1 and conversely, N15N14+N15(t=0)=0

(4) the depletion ratios vanish in finite time, i.e., N14N14+N15(t→∞)=0 and conversely, N15N14+N15(t→∞)=1

(5) the (m + 0) and (m + 1) isotopologues are assumed to be measures for the abundance of N14 and N15 isotopes respectively.

The model is thus represented as follows:

Given,d[N15]dt=kcat+∗([N]total−[N15])

Upon integration,ln([N]total[N]total−[N15])=kcat+×t

N15(t)=[N]total∗(1−e−kcat+t)=(N14+N15)×(1−e−kcat+t)

Hence,fraction(t)=N15N14+N15(t)=1−e−kcat+t

This approach assumes a direct correlation between the rate of 15N incorporation and the metabolic turnover rates, allowing us to infer the in vivo catalytic speed and efficiency of sphingolipid biosynthesis reactions. The linear regression of the transformed nonlinear model curve facilitates the quantification of kcat+ or TOF values, providing insights into the metabolic demands and efficiency of sphingolipid biosynthesis in the cultures. However, we used a widely known technique called bootstrapping to handle data variability, which allows us to estimate the uncertainties.76 We do this by assuming that the variations we see in our experiments follow a normal (Gaussian) distribution pattern. This means we’re deliberately adding a layer of uncertainty to better reflect the natural fluctuations seen in our measurements (per culture C and D), especially regarding the ratios’ bounds we’ve calculated. As to the ratios’ bounds, as shown in Figure 2C, we used four technical replicates for each labeling condition: 4 replicates from cultures C and D for 15N labeling, and 4 replicates from cultures A and B for 14N labeling. This includes two technical replicates per biological duplicate (2 for cultures C, 2 for culture D, 2 for culture A, 2 for culture B) to obtain the mean and bounds for the enrichment profiles.

Estimating dynamic flux using regularized DMFA

A regularized form of the dynamic metabolic flux analysis (L-DMFA) framework by Leighty and Antoniewicz (2011)77 was used to predict the dynamic metabolic fluxes, v(t). DMFA methodology uses a DOA approach to fit concentration measurements directly and relies on non-steady state mass balances for metabolite pools. The model non-steady state mass balances for metabolite pools are represented by:(Equation 1) dcdt=S·K·u(t)with,K=null(Sbal)

In this context and for a given metabolic network, we denote, S as the stoichiometric matrix, K as the null space matrix of S, and c(t) as the vector of metabolite concentrations. Metabolites are categorized as balanced, typically intracellular, or non-balanced, typically external. Leighty and Antoniewicz, (2011) defined and determined internal free fluxes, u(t) as a function of linear splines. It divides the time domain [t1,tn] of a culture into smaller intervals [tk−1,tk], assuming constant rate changes in fluxes between distinct DMFA time points. Hence, the flux dynamics are modeled as linear flux changes between two points as:(Equation 2) v(t)=K·uk−1·(1−t−tk−1tk−tk−1)+K·uk·(t−tk−1tk−tk−1)

In a more compact form:(Equation 3) v(t)=K·U·κ(t,tk)

where, U=[uk−1,uk,…,un], uk−1 and uk denote the free internal fluxes to be fitted at times tk−1 and tk respectively. κ(t,tk) is a time dependent parameter matrix for the linear spline with its collocations defined in.77 Integrating Equations 1 and 3 yields:(Equation 4) c(t)=c0+S·K·U·γ(t,tk)

where, c0 is a fitted variable defined as the initial metabolite concentrations, and γ(t,tk) is the integral form of κ(t,tk). For estimating flux transients, the variance-weighted sum of squared residuals (SSR) is minimized subject to Equation 4 as follows:(Equation 5) minSSR=∑(ci−ci,m)2σi2Subjectto:ci=ci,0+Si·K·U·γ,∀i∈I

where, ci,m are measured concentration profiles for metabolite pools; I are sets of metabolites in the model; σi are the specific measurement variances assumed as experimentally determined fluctuations. Hence, direct linear explicit and global optimums, p are obtained from:(Equation 6) p∶=[c0U]=H−1·J

where, the Hessian and Jacobian matrices, H and J respectively are calculated as:H≝∑[(δciδp)T·σi2·(δciδp)]andJ≝∑[(δciδp)T·σi2·ci,m]

with the derivatives of ci,δciδp≝[Ei(γT·Di)]

where, Di=diag(Si·K) and Ei∋δi,j={1ifi=j0ifi≠j

With the current formulation, solving the linear system involves inverting the Hessian matrix, H to estimating the parameter set, p. However, when data includes noisy (such as in our study the large variation in the distributions of σi) or incomplete measurements the problem might be ill-posed; that is, it may not have a unique solution, or small changes in the data might lead to drastically different solutions. H could as well be near singular or ill-conditioned, its inversion can lead to numerical instability. Overfitting is another concern, where the estimated model fits the noise in the data rather than the underlying process. Tikhonov regularization of H was used to address these issues by adding a regularization term that conditions H, making its inversion more stable. This regularization promotes numerical stability and generalization The regularized Hessian, H′ is given by:(Equation 7) H′=H+λI

where, I is an identity matrix. With the regularized Hessian, H′, the regularization parameter, λ, balances the trade-off between fitting the data and regularizing the solution in Equation 6, and can be chosen via methods like cross-validation or information criteria. In this study, the regularization parameter, λ, is chosen through a systematic grid search optimization utilizing the SSR, and Bayesian Information Criterion (BIC). Optimal λ was chosen with the smallest BIC value. In our study, we obtained λopt=1×10−4 for both independent biological duplicates.

Building upon the L-DMFA method, which assumes linear flux changes between designated DMFA time points, our approach in the same fashion employs an iterative process to determine the optimal set of linear time intervals. This process aims to minimize the Sum of Squared Residuals (SSR), (see Figure S3) thereby ensuring the model’s convergence and bolstering its predictive accuracy. However, central, and specifically to our approach is enhancing the L-DMFA fitting process by the use of cubic spline interpolation. This technique not only incorporates every data measurement point but also generates additional points from cubic splines within these intervals, with the stipulation that each interval encompasses strictly one data point, no less, no more. These supplementary points serve as inflection points, enabling the modeling of linear flux changes between them. The result is a significantly smoother curve that more faithfully represents the experimental data’s non-linear flux changes, thus improving the model’s precision. This is coupled with earlier regularization steps for metabolic flux analysis. Upon this refinement, the formulation aligns with the L-DMFA approach as delineated by Leighty and Antoniewicz (2011).77 The predicted fluxes, based on the average concentration profiles, are intended to reflect the typical flux distributions for cultures C and D. Consequently, we have roughly estimated the standard error for these average flux distributions, as illustrated in Figure 5. Following this, we adopted a comprehensive strategy that involves uniform convex polytope sampling, referred to as dynamic flux sampling. This technique was employed to achieve a more accurate distribution of fluxes, which is elaborated upon in the subsequent two subsections.

Additionally, we deliberately maintained the reversibility of all reactions in this pathway because we are analyzing data from inside cells. This decision was made to expand the range of our predictions, allowing for a wider exploration of metabolic pathways without being limited to predetermined outcomes. The regularized DMFA framework, r-DMFA as a modified version of the MATLAB implementation developed by Leighty and Antoniewicz (2011)77 was scripted and used to calculate the reaction rates for each enzymatic reaction in the network.

Estimating enzyme metabolic demands

By intertwining thermodynamics and enzyme kinetics, we used the concept of enzyme metabolic demands which facilitates a holistic understanding of underlying mechanism dictating the efficiency, speed, and equilibrium of enzymatic reactions in our network. This proffers an estimation of active enzyme molecules partaking in the reaction. To this end, we determined the variable enzyme metabolic demands, EMD for each reaction in the network as defined by Noor et al.,.54 This was summed over all reactions to also obtain the enzyme level metabolic demand. For a single reaction, j, the enzyme metabolic demand is defined as the ratio of metabolic reaction flux, vj(t) and the relative/apparent turnover frequency (kcat+), our experimentally fitted TOF and a measure of reaction speed, given by the equation:(Equation 8) EMD,j(t)=vj(t)kcat+where,kcat+=kcat+∗·γ·Κ

kcat+∗ will be the actual enzyme turnover number independent of thermodynamic and substrate saturation. γ∈[−1,1], is the thermodynamic driving force. Κ∈[0,1], is the enzyme specific substrate saturation factor.

Analyzing transient flux distribution through dynamic flux sampling

A concerted approach employing a uniform convex polytope sampling strategy was adopted to delve into the high-dimensional dynamic flux solution spaces, encompassing diverse alternate optimal dynamic flux distributions. This entailed uniformly generating random transient concentration profiles within the known experimental constraints, leading to the truncation of the bounded flux polytope—a modality prevalently utilized for characterizing metabolic solution spaces.57,58 To enable an efficient sampling process and mitigate the inefficiencies inherent in rejection sampling, the r-DMFA framework was further modified by incorporating a global optimization strategy facilitated by the BARON global solver.59

Following this, a flux modulation ratio was computed for each culture, defined as the fold changes of each flux in the generated ensemble dynamic model relative to the average representative flux distribution per culture. These flux fold changes served to characterize the influence of each reaction on the variations observed within the ensemble of dynamic flux modulation ratios, providing insight into the reaction-driven flux alterations across different cultures.

In the characterization of the emergent dynamic flux fold changes, we employed the two-dimensional t-Distributed Stochastic Neighbor Embedding (t-SNE) for dimensionality reduction. t-SNE, a nonlinear technique, excels in the visualization of high-dimensional datasets within two or three-dimensional spaces.78 This method transforms the Euclidean distances in high-dimensional space into conditional probabilities that mirror similarities among data points, thereby facilitating an intuitive representation of complex data structures. Subsequent to the application of t-SNE, the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm was utilized to delineate representative clusters within the resultant two-dimensional space.79 This clustering technique enables the identification of dense regions of data points, which are indicative of relative clusters, thereby augmenting the analysis with a robust spatial interpretation of the data’s inherent groupings.

To elucidate the relationships and correlations between the predicted metabolic concentration profiles and those obtained through dynamic sampling, Spearman’s rank correlation analysis was employed. This statistical approach, leveraging the robustness of Spearman’s non-parametric method, allowed for the assessment of monotonic relationships between the two sets of data. By comparing the ranked positions of metabolic concentrations across both predicted and dynamically sampled profiles, we were able to quantify the degree of correlation without the assumption of linear relationships or normal data distribution. MATLAB’s “corr’ function facilitated this analysis, providing not only the Spearman’s rank correlation coefficients but also the corresponding p-values. These p-values, derived either from exact permutation distributions for smaller datasets or from large-sample approximations for larger ones, offered a statistical basis to gauge the significance of the correlations observed. This approach enabled a rigorous evaluation of the predictive accuracy of our dynamic models against the backdrop of experimentally derived data, underscoring the potential variances and consistencies in metabolic behavior.

The distribution of Spearman’s correlation coefficients, rho (ρ), and their respective p-values were meticulously analyzed to discern sensitivity of approach and statistically significant correlations, thereby illuminating the strength and direction of the relationships between predicted and sampled metabolic concentration profiles. This statistical scrutiny was paramount in distinguishing genuine biological correlations from those potentially arising by chance. The insights garnered from this analysis not only augmented our understanding of the metabolic flux distributions but also provided a foundation for further exploration into the dynamic intricacies of metabolic pathways. Through this methodical examination, we aimed to broaden the scope of our investigation, facilitating a more comprehensive understanding of the metabolic dynamics at play within various cultures and even under varying experimental conditions, which could be the case with high data variability.

Analyzing the sensitivity and robustness of predictions

In analyzing the sensitivity and robustness of predictions within the context of metabolic pathway regulation, as depicted in Figure S6 which offers significant insights through the analysis of flux modulation ratios. These ratios, defined as fold changes of flux in relation to the average flux per culture, have been computed for each dynamic flux ensemble and culture, effectively characterizing the influence and variability of individual reactions. The application of t-Distributed Stochastic Neighbor Embedding (t-SNE) for dimensionality reduction has successfully transformed high-dimensional flux data into a two-dimensional space, where the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm then identified distinct clusters indicative of reaction partitioning within the metabolic network.

The sensitivity and robustness of the predictions enabled by DBSCAN clustering process across different sample sizes in the analysis of metabolic flux distributions is shown in Panel A of Figure S5. The panel displays the DBSCAN clusters for various reactions in both C and D cultures, with sample sizes ranging from 1E+4 to 4E+6 and a separate category for outliers. Outliers here are determined by the Spearman’s rho (ρ) correlation distributions with respect to the fits of the concentration profiles as would be in Figure S4. Note that this follows the definition that a value that is more than three scaled median absolute deviations (MAD) from the median of the Spearman’s rho (ρ) correlation values. Moreover, the color-coding represents different clusters, while the numbers indicate the cluster to which each reaction belongs or whether it is classified to have been impacted by noise effects (i.e., does not belong to any distinct cluster).

The consistent clustering of reactions across varying sample sizes demonstrates the robustness of the predictive model and approach. Reactions that maintain their cluster assignment despite changes in sample size indicate a strong and stable signal in the predictions, which could reflect the reactions’ consistent behavior under different conditions. For instance, reactions that consistently fall within the red cluster across different sample sizes suggest a high level of regulation and a critical role in maintaining cellular homeostasis, as these reactions are less sensitive to fluctuations in sample size, indicating a stable role in metabolic processes.

Conversely, reactions that shift between clusters or to the outlier category as the sample size changes could either indicate a sensitivity to the data density or reveal insights into the biological variability of these reactions. This observation underlines the capacity of the model to differentiate between tightly regulated and more flexible metabolic reactions, thereby highlighting the nuanced balance within the metabolic network that sustains cellular function under various conditions.

The minimal variation in cluster assignments, especially for significant reactions across the full range of sample sizes, underscores the method’s precision. This consistent clustering across different sample sizes also supports the robustness of the metabolic model, suggesting that the predictive insights it offers are reliable and can withstand the test of varying data quantities, an essential feature and focus of this study.

Assessing predicted metabolic concentration profiles

Residuals analysis

We assessed the accuracy of our predicted dynamic profiles by analyzing the residuals, which demonstrated a tendency to cluster around zero, indicating a good fit for the species within both C and D cultures (refer to Figure S3 for details).

Correlation analysis

To evaluate the relationship between predicted metabolic concentration profiles and dynamically sampled profiles, we employed Spearman’s rank correlation coefficient. This non-parametric measure assesses how well the relationship between two variables can be described using a monotonic function. Unlike Pearson’s correlation, Spearman’s correlation does not assume a linear relationship between variables or that the variables are normally distributed. Instead, it relies on the ranked values of the data points, making it more robust against non-normal data distributions and outliers. This methodology facilitated a comprehensive and statistically robust analysis of the correlations between predicted time metabolic concentration profiles and the randomly sampled dynamic profiles, broadening the scope of our investigation into metabolic pathway dynamics.

Spearman’s rank correlation coefficient - mathematical overview

Spearman’s rank correlation coefficient, rho (ρ), is defined as:ρ=1−6∑dj2n(n2−1)

where:

dj is the difference between the ranks of corresponding variables cipred(j) and cisampled(j) for j∈[t1,tn]

n is the number of observations, i.e., time points, j.

A rho (ρ) value of +1 indicates a perfect positive monotonic relationship, −1 a perfect negative monotonic relationship, and 0 indicates no monotonic relationship.

Implementation and Statistical Significance

We utilized MATLAB’s `corr` function to compute Spearman’s rho (ρ) and associated p-values between the predicted and sampled metabolic concentration profiles. The `corr` function also offers the capability to calculate Kendall’s tau, but for our analysis, Spearman’s rho (ρ) was chosen for its ability to provide a measure of monotonic relationships. Depending on the size of our datasets, the function automatically selects between exact permutation distributions (for small sample sizes) and large-sample approximations to compute p-values. These p-values are crucial for determining the statistical significance of the correlations we observed. A lower p-value, particularly one less than 0.05, indicates that the correlation (ρ) between two variables significantly deviates from zero, suggesting a meaningful association contrary to the null hypothesis, which posits no relationship.

Generating perturbation models by knockouts (KO)

Perturbations were induced through single enzyme and reaction knockouts, performed by leveraging additional constraints to the r-DMFA framework, grounded on the Minimization of Metabolic Adjustment (MOMA) principle,60 to predict dynamic model responses. In other words, our MOMA-based r-DMFA framework employed quadratic programming, relaxing optimal WT flux assumptions, facilitating precise metabolic phenotype predictions of mutants, and elucidating cellular adaptation to enzymatic dysfunction. Furthermore, this approach utilized convex transformation methods in the global optimization algorithm of the BARON solver.61 Moreover, the perturbations could take on the form of enzyme overexpression, inhibition, or knockout and is formulated as an additional constraint to the problem in Equation 5 given by:(Equation 9) ei,1˜vj,min(t)≤vj(t)≤vj,max(t)ei,1˜

where, ei,1˜ represents the introduced regulation fold change in enzyme expression compared to the reference state.

After each perturbation, the goal was to compare the predicted responses from single knockouts against the experimental WT predictions. An evaluation metric, zd, and spearman’s rank correlation was employed to evaluate the metabolic distance between each mutant dynamic flux phenotype and the WT dynamic flux phenotype. A supplementary module was integrated into the original framework to facilitate these modifications. Additionally, a modified version of the MATLAB implementation as developed by Leighty and Antoniewicz (2011)77 was utilized to conduct the simulations in this study. The dynamic model along with all requisite scripts for generating the results of this study are made accessible via the link: https://github.com/ssbio/r-DMFA.

Quantifying metabolic distances in dynamic perturbation models

An evaluation metric, defined as the average relative deviation between WT determined fluxes and those predicted by perturbed (mutant) dynamic models, was employed. The metric utilizes the coefficient of variation80 to scale the terms in the metric function, zd, thereby capturing potential degenerate polytopic uncertainty associated with experimental measurements. Consequently, reactions with a tighter confidence interval contribute more significantly to the function. The metric is expressed as:(Equation 10) zd=∑j∈Nr∑t∈Tln(1CVj(t)|vj(t)−vjWT(t)vjWT(t)|)

where, vjWT(t) represents the WT flux profile of reaction j; CVj(t) denotes the coefficient of variation for each mutant dynamic flux of reaction j at each timepoint, employed to assess the relative dispersion of the flux-space solution with sensitivity to experimental measurement uncertainty. NT represents the number of considered timepoints of T, the time domain [t1,…,tn] ; Nr is the set of nonzero flux reactions in the WT dynamic model. It is noteworthy that reactions that do not carry flux for the WT may exhibit nonzero flux in a perturbed dynamic model. Excluding zero-flux WT reactions from the comparative measure of metabolic distance enables a rigorous and relative comparison across mutants. Mathematically, this excludes jump discontinuities.

Additionally, to assess the metabolic differences between mutant and WT strains, Spearman’s rank correlation was utilized. This method calculates Spearman’s rank correlation coefficients for each metabolic reaction across various time points. Accompanying each output correlation coefficient, designated as rho, is a statistical measure, the p-value, indicating statistical significance. From this, two primary metrics were derived: the ratios of positive correlation reactions and negative correlation reactions. These metrics represent the proportion of reactions displaying significant positive (rho >0.5, p < 0.05) and negative (rho < −0.5, p < 0.05) correlations, respectively. This provides a succinct depiction of metabolic congruence or divergence between the two strains.

Quantification and statistical analysis

Within this study, various statistical and analytical tools were employed. The two-sample Kolmogorov-Smirnov test (KS test) (Frank J. Massey, 1951)56 was utilized to compare distributions, while the Bhattacharya coefficient was employed for overlap measurement between distributions. Bootstrap resampling was employed to provide estimation accuracy, and linear regression was used for examining relationships between variables. Additionally, Spearman’s rank correlation coefficient was implemented to evaluate the monotonic relationships between dynamic flux distributions, particularly given that the data in this study did not satisfy the requirements of linearity and homoscedasticity, thereby providing a non-parametric alternative to assess associations without assuming a specific distribution.

Supplemental information

Document S1. Figures S1–S6, Tables S1, and S2

Data S1. Temporal abundance of 15N in sphingolipid species, related to Figure 2D

Acknowledgments

A.O. wishes to express gratitude for the financial support provided for this research. This includes the 10.13039/100000002 National Institutes of Health (NIH) MIRA Award (5R35GM143009 ) awarded to R.S. and the 10.13039/100000001 National Science Foundation (NSF) grant (MCB 1818297 ) awarded to E.B.C. and R.S. We gratefully acknowledge the source code to the initial L-DMFA framework provided by Maciek Antoniewicz from 10.13039/100007270 University of Michigan .

Author contributions

Conceptualization, R.S. and E.B.C.; methodology, A.O.; software, A.O.; investigation, A.O., A.G.S., and R.E.C.; resources, R.S. and E.B.C.; data curation, A.O.; writing—original draft, A.O.; writing—review and editing, A.G.S., R.E.C., A.A., E.B.C., and R.S.; visualization, A.O.; supervision, R.S. and E.B.C.; funding acquisition, E.B.C. and R.S.; project administration, E.B.C. and R.S.

Declaration of interests

The authors declare no competing interests.

Supplemental information can be found online at https://doi.org/10.1016/j.isci.2024.110675.
==== Refs
References

1 Snider J.M. Snider A.J. Obeid L.M. Luberto C. Hannun Y.A. Probing de novo sphingolipid metabolism in mammalian cells utilizing mass spectrometry J. Lipid Res. 59 2018 1046 1057 10.1194/JLR.D081646 29610123
2 Michaelson L.V. Napier J.A. Molino D. Faure J.D. Plant sphingolipids: Their importance in cellular organization and adaption Biochim. Biophys. Acta 1861 2016 1329 1335 10.1016/J.BBALIP.2016.04.003 27086144
3 Buré C. Cacas J.L. Mongrand S. Schmitter J.M. Characterization of glycosyl inositol phosphoryl ceramides from plants and fungi by mass spectrometry Anal. Bioanal. Chem. 406 2014 995 1010 10.1007/S00216-013-7130-8 23887274
4 Luttgeharm K.D. Kimberlin A.N. Cahoon E.B. Plant sphingolipid metabolism and function Subcell. Biochem. 86 2016 249 286 10.1007/978-3-319-25979-6_11 27023239
5 Worrall D. Liang Y.K. Alvarez S. Holroyd G.H. Spiegel S. Panagopulos M. Gray J.E. Hetherington A.M. Involvement of sphingosine kinase in plant cell signalling Plant J. 56 2008 64 72 10.1111/J.1365-313X.2008.03579.X 18557834
6 Magnin-Robert M. Le Bourse D. Markham J. Dorey S. Clément C. Baillieul F. Dhondt-Cordelier S. Modifications of Sphingolipid Content Affect Tolerance to Hemibiotrophic and Necrotrophic Pathogens by Modulating Plant Defense Responses in Arabidopsis Plant Physiol. 169 2015 2255 2274 10.1104/PP.15.01126 26378098
7 Merrill A.H. De Novo Sphingolipid Biosynthesis: A Necessary, but Dangerous, Pathway J. Biol. Chem. 277 2002 25843 25846 10.1074/JBC.R200009200 12011104
8 Merrill A.H. Sphingolipid and glycosphingolipid metabolic pathways in the era of sphingolipidomics Chem. Rev. 111 2011 6387 6422 10.1021/CR2002917 21942574
9 Hanada K. Serine palmitoyltransferase, a key enzyme of sphingolipid metabolism Biochim. Biophys. Acta 1632 2003 16 30 10.1016/S1388-1981(03)00059-3 12782147
10 Markham J.E. Lynch D.V. Napier J.A. Dunn T.M. Cahoon E.B. Plant sphingolipids: function follows form Curr. Opin. Plant Biol. 16 2013 350 357 10.1016/J.PBI.2013.02.009 23499054
11 Markham J.E. Molino D. Gissot L. Bellec Y. Hématy K. Marion J. Belcram K. Palauqui J.C. Satiat-Jeunemaître B. Faure J.D. Sphingolipids Containing Very-Long-Chain Fatty Acids Define a Secretory Pathway for Specific Polar Plasma Membrane Protein Targeting in Arabidopsis Plant Cell 23 2011 2362 2378 10.1105/TPC.110.080473 21666002
12 Ternes P. Feussner K. Werner S. Lerche J. Iven T. Heilmann I. Riezman H. Feussner I. Disruption of the ceramide synthase LOH1 causes spontaneous cell death in Arabidopsis thaliana New Phytol. 192 2011 841 854 10.1111/j.1469-8137.2011.03852.x 21883234
13 Msanne J. Chen M. Luttgeharm K.D. Bradley A.M. Mays E.S. Paper J.M. Boyle D.L. Cahoon R.E. Schrick K. Cahoon E.B. Glucosylceramides are critical for cell-type differentiation and organogenesis, but not for cell viability in Arabidopsis Plant J. 84 2015 188 201 10.1111/TPJ.13000 26313010
14 Greenberg A.J. Hackett S.R. Harshman L.G. Clark A.G. Environmental and genetic perturbations reveal different networks of metabolic regulation Mol. Syst. Biol. 7 2011 563 10.1038/MSB.2011.96 22186737
15 Nair A. Chauhan P. Saha B. Kubatzky K.F. Conceptual Evolution of Cell Signaling Int. J. Mol. Sci. 20 2019 3292 10.3390/IJMS20133292
16 Guo L. Wang X. Crosstalk between phospholipase D and sphingosine kinase in plant stress signaling Front. Plant Sci. 3 2012 51 10.3389/FPLS.2012.00051
17 Qin X. Zhang R.X. Ge S. Zhou T. Liang Y.K. Sphingosine kinase AtSPHK1 functions in fumonisin B1-triggered cell death in Arabidopsis Plant Physiol. Biochem. 119 2017 70 80 10.1016/J.PLAPHY.2017.08.008 28846870
18 Imai H. Nishiura H. Phosphorylation of Sphingoid Long-chain Bases in Arabidopsis: Functional Characterization and Expression of the First Sphingoid Long-chain Base Kinase Gene in Plants Plant Cell Physiol. 46 2005 375 380 10.1093/pcp/pci023 15695468
19 Shi L. Bielawski J. Mu J. Dong H. Teng C. Zhang J. Yang X. Tomishige N. Hanada K. Hannun Y.A. Zuo J. Involvement of sphingoid bases in mediating reactive oxygen intermediate production and programmed cell death in Arabidopsis Cell Res. 17 2007 1030 1040 10.1038/cr.2007.100 18059378
20 Ray P.D. Huang B.W. Tsuji Y. Reactive oxygen species (ROS) homeostasis and redox regulation in cellular signaling Cell. Signal. 24 2012 981 990 10.1016/J.CELLSIG.2012.01.008 22286106
21 Aguilera-Romero A. Gehin C. Riezman H. Sphingolipid homeostasis in the web of metabolic routes Biochim. Biophys. Acta 1841 2014 647 656 10.1016/j.bbalip.2013.10.014 24184514
22 Markham J.E. Jaworski J.G. Rapid measurement of sphingolipids from Arabidopsis thaliana by reversed-phase high-performance liquid chromatography coupled to electrospray ionization tandem mass spectrometry Rapid Commun. Mass Spectrom. 21 2007 1304 1314 10.1002/RCM.2962 17340572
23 Haynes C.A. Allegood J.C. Wang E.W. Kelly S.L. Sullards M.C. Merrill A.H. Factors to consider in using [U-13C]palmitate for analysis of sphingolipid biosynthesis by tandem mass spectrometry J. Lipid Res. 52 2011 1583 1594 10.1194/jlr.D015586 21586681
24 Wigger D. Gulbins E. Kleuser B. Schumacher F. Monitoring the Sphingolipid de novo Synthesis by Stable-Isotope Labeling and Liquid Chromatography-Mass Spectrometry Front. Cell Dev. Biol. 7 2019 210 10.3389/FCELL.2019.00210
25 Siow D. Sunkara M. Morris A. Wattenberg B. Regulation of de novo sphingolipid biosynthesis by the ORMDL proteins and sphingosine kinase-1 Adv. Biol. Regul. 57 2015 42 54 10.1016/J.JBIOR.2014.09.002 25319495
26 You L. Zhang B. Tang Y.J. Application of Stable Isotope-Assisted Metabolomics for Cell Metabolism Studies Metabolites 4 2014 142 165 10.3390/METABO4020142 24957020
27 Chen X. Alonso A.P. Shachar-Hill Y. Dynamic metabolic flux analysis of plant cell wall synthesis Metab. Eng. 18 2013 78 85 10.1016/J.YMBEN.2013.04.006 23644173
28 Ecker J. Liebisch G. Application of stable isotopes to investigate the metabolism of fatty acids, glycerophospholipid and sphingolipid species Prog. Lipid Res. 54 2014 14 31 10.1016/J.PLIPRES.2014.01.002 24462586
29 Allen D.K. Young J.D. Tracing metabolic flux through time and space with isotope labeling experiments Curr. Opin. Biotechnol. 64 2020 92 100 10.1016/J.COPBIO.2019.11.003 31864070
30 Martínez-Montañés F. Schneiter R. Tools for the analysis of metabolic flux through the sphingolipid pathway Biochimie 130 2016 76 80 10.1016/J.BIOCHI.2016.05.009 27208414
31 Kuchař L. Asfaw B. Rybová J. Ledvinová J. Tandem Mass Spectrometry of Sphingolipids: Applications for Diagnosis of Sphingolipidoses Adv. Clin. Chem. 77 2016 177 219 10.1016/BS.ACC.2016.06.004 27717417
32 Li S. Park Y. Duraisingham S. Strobel F.H. Khan N. Soltow Q.A. Jones D.P. Pulendran B. Predicting Network Activity from High Throughput Metabolomics PLoS Comput. Biol. 9 2013 e1003123 10.1371/JOURNAL.PCBI.1003123
33 Allen D.K. Bates P.D. Tjellström H. Tracking the metabolic pulse of plant lipid production with isotopic labeling and flux analyses: Past, present and future Prog. Lipid Res. 58 2015 97 120 10.1016/J.PLIPRES.2015.02.002 25773881
34 Schroeder W.L. Saha R. Introducing an Optimization- and explicit Runge-Kutta- based Approach to Perform Dynamic Flux Balance Analysis Sci. Rep. 10 2020 9241 10.1038/S41598-020-65457-4
35 Antoniewicz M.R. 13C metabolic flux analysis: optimal design of isotopic labeling experiments Curr. Opin. Biotechnol. 24 2013 1116 1121 10.1016/J.COPBIO.2013.02.003 23453397
36 Ma F. Jazmin L.J. Young J.D. Allen D.K. Isotopically nonstationary 13C flux analysis of changes in Arabidopsis thaliana leaf metabolism due to high light acclimation Proc. Natl. Acad. Sci. USA 111 2014 16967 16972 10.1073/PNAS.1319485111 25368168
37 Wall M.L. Pound L.D. Trenary I. O’Brien R.M. Young J.D. Novel stable isotope analyses demonstrate significant rates of glucose cycling in mouse pancreatic islets Diabetes 64 2015 2129 2137 10.2337/DB14-0745 25552595
38 Li L. Nelson C.J. Trösch J. Castleden I. Huang S. Millar A.H. Protein degradation rate in Arabidopsis thaliana leaf growth and development Plant Cell 29 2017 207 228 10.1105/TPC.16.00768 28138016
39 Antoniewicz M.R. Dynamic metabolic flux analysis — tools for probing transient states of metabolic networks Curr. Opin. Biotechnol. 24 2013 973 978 10.1016/J.COPBIO.2013.03.018 23611566
40 Hackett S.R. Zanotelli V.R.T. Xu W. Goya J. Park J.O. Perlman D.H. Gibney P.A. Botstein D. Storey J.D. Rabinowitz J.D. Systems-level analysis of mechanisms regulating yeast metabolic flux Science (1979) 354 2016 aaf2786 10.1126/SCIENCE.AAF2786
41 Di Filippo M. Pescini D. Galuzzi B.G. Bonanomi M. Gaglio D. Mangano E. Consolandi C. Alberghina L. Vanoni M. Damiani C. INTEGRATE: Model-based multi-omics data integration to characterize multi-level metabolic regulation PLoS Comput. Biol. 18 2022 e1009337 10.1371/JOURNAL.PCBI.1009337
42 Johnson C.H. Ivanisevic J. Siuzdak G. Metabolomics: beyond biomarkers and towards mechanisms Nat. Rev. Mol. Cell Biol. 17 2016 451 459 10.1038/nrm.2016.25 26979502
43 Cahoon R.E. Solis A.G. Markham J.E. Cahoon E.B. Mass spectrometry-based profiling of plant sphingolipids from typical and aberrant metabolism Methods Mol. Biol. 2295 2021 157 177 10.1007/978-1-0716-1362-7_10 34047977
44 Kasumov T. Huang H. Chung Y.M. Zhang R. McCullough A.J. Kirwan J.P. Quantification of ceramide species in biological samples by liquid chromatography electrospray ionization tandem mass spectrometry Anal. Biochem. 401 2010 154 161 10.1016/J.AB.2010.02.023 20178771
45 Yamaguchi M. Miyashita Y. Kumagai Y. Kojo S. Change in liver and plasma ceramides during D D-galactosamine-induced acute hepatic injury by LC-MS/MS Bioorg. Med. Chem. Lett. 14 2004 4061 4064 10.1016/j.bmcl.2004.05.046 15225726
46 Baranowski M. Blachnio A. Zabielski P. Gorski J. Pioglitazone induces de novo ceramide synthesis in the rat heart Prostaglandins Other Lipid Mediat. 83 2007 99 111 10.1016/j.prostaglandins.2006.10.004 17259076
47 Scherer M. LeuthÃ K. Ecker J. Schmitz G. Liebisch G. A rapid and quantitative LC-MS/MS method to profile sphingolipids 10.1194/jlr.D005322 2010
48 Hebing L. Neymann T. Engell S. Application of dynamic metabolic flux analysis for process modeling: Robust flux estimation with regularization, confidence bounds, and selection of elementary modes Biotechnol. Bioeng. 117 2020 2058 2073 10.1002/BIT.27340 32196640
49 Alsiyabi A. Solis A.G. Cahoon E.B. Saha R. Dissecting the regulatory roles of ORM proteins in the sphingolipid pathway of plants PLoS Comput. Biol. 17 2021 e1008284 10.1371/JOURNAL.PCBI.1008284
50 Bhattacharyya A. On a measure of divergence between two multinomial populations Sankhya 7 1946 401 406
51 Bhattacharyya A. On a Measure of Divergence Between Two Statistical Populations Defined by their Probability Distributions Bull. Calcutta Math. Soc. 35 1943 99 110
52 Liu N.-J. Hou L.-P. Bao J.-J. Wang L.-J. Chen X.-Y. N-j L. L-p H. J-j B. L-j W. X-y C. Sphingolipid metabolism, transport, and functions in plants: Recent progress and future perspectives Plant Commun. 2 2021 100214 10.1016/j.xplc.2021.100214
53 Markham J.E. Li J. Cahoon E.B. Jaworski J.G. Separation and Identification of Major Plant Sphingolipid Classes from Leaves J. Biol. Chem. 281 2006 22684 22694 10.1074/JBC.M604050200 16772288
54 Noor E. Flamholz A. Bar-Even A. Davidi D. Milo R. Liebermeister W. The Protein Cost of Metabolic Fluxes: Prediction from Enzymatic Rate Laws and Cost Minimization PLoS Comput. Biol. 12 2016 e1005167 10.1371/JOURNAL.PCBI.1005167
55 Noor E. Flamholz A. Liebermeister W. Bar-Even A. Milo R. A note on the kinetics of enzyme action: A decomposition that highlights thermodynamic effects FEBS Lett. 587 2013 2772 2777 10.1016/j.febslet.2013.07.028 23892083
56 Massey F.J. The Kolmogorov-Smirnov Test for Goodness of Fit J. Am. Stat. Assoc. 46 1951 68 78 10.1080/01621459.1951.10500769
57 Herrmann H.A. Dyson B.C. Vass L. Johnson G.N. Schwartz J.M. Flux sampling is a powerful tool to study metabolism under changing environmental conditions NPJ Syst. Biol. Appl. 2019 32 38 10.1038/s41540-019-0109-0 31482008
58 Schellenberger J. Palsson B.Ø. Use of Randomized Sampling for Analysis of Metabolic Networks J. Biol. Chem. 284 2009 5457 5461 10.1074/JBC.R800048200 18940807
59 Amaran S. Sahinidis N.V. Global optimization of nonlinear least-squares problems by branch-and-bound and optimality constraints TOP 20 2012 154 172 10.1007/S11750-011-0178-8
60 Segrè D. Vitkup D. Church G.M. Analysis of optimality in natural and perturbed metabolic networks Proc. Natl. Acad. Sci. USA 99 2002 15112 15117 10.1073/PNAS.232349399 12415116
61 Nohra C.J. Sahinidis N.V. Global optimization of nonconvex problems with convex-transformable intermediates J. Glob. Optim. 72 2018 255 276 10.1007/S10898-018-0631-4
62 Schellenberger J. Zielinski D.C. Choi W. Madireddi S. Portnoy V. Scott D.A. Reed J.L. Osterman A.L. Palsson B.T. Predicting outcomes of steady-state 13C isotope tracing experiments using Monte Carlo sampling BMC Syst. Biol. 6 2012 1 14 10.1186/1752-0509-6-9 22222070
63 Chen M. Markham J.E. Dietrich C.R. Jaworski J.G. Cahoon E.B. Sphingolipid Long-Chain Base Hydroxylation Is Important for Growth and Regulation of Sphingolipid Content and Composition in Arabidopsis Plant Cell 20 2008 1862 1878 10.1105/TPC.107.057851 18612100
64 Pralhada Rao R. Vaidyanathan N. Rengasamy M. Mammen Oommen A. Somaiya N. Jagannath M.R. Sphingolipid Metabolic Pathway: An Overview of Major Roles Played in Human Diseases J. Lipids 2013 2013 1 12 10.1155/2013/178910
65 Hannun Y.A. Obeid L.M. Principles of bioactive lipid signalling: lessons from sphingolipids Nat. Rev. Mol. Cell Biol. 9 2008 139 150 10.1038/nrm2329 18216770
66 González Solís A. Dissecting the Regulatory Network of Sphingolipid Biosynthesis in Plants PhD Dissertation 2020 ETD collection for University of Nebraska-Lincoln AAI28259034 https://digitalcommons.unl.edu/dissertations/AAI28259034
67 Liang H. Yao N. Song J.T. Luo S. Lu H. Greenberg J.T. Ceramides modulate programmed cell death in plants Genes Dev. 17 2003 2636 2641 10.1101/GAD.1140503 14563678
68 Ukawa T. Banno F. Ishikawa T. Kasahara K. Nishina Y. Inoue R. Tsujii K. Yamaguchi M. Takahashi T. Fukao Y. Sphingolipids with 2-hydroxy fatty acids aid in plasma membrane nanodomain organization and oxidative burst Plant Physiol. 189 2022 839 857 10.1093/PLPHYS/KIAC134 35312013
69 Chen M. Markham J.E. Cahoon E.B. Sphingolipid Δ8 unsaturation is important for glucosylceramide biosynthesis and low-temperature performance in Arabidopsis Plant J. 69 2012 769 781 10.1111/J.1365-313X.2011.04829.X 22023480
70 Nishikawa M. Hosokawa K. Ishiguro M. Minamioka H. Tamura K. Hara-Nishimura I. Takahashi Y. Shimazaki K.I. Imai H. Degradation of Sphingoid Long-Chain Base 1-Phosphates (LCB-1Ps): Functional Characterization and Expression of AtDPL1 Encoding LCB-1P Lyase Involved in the Dehydration Stress Response in Arabidopsis Plant Cell Physiol. 49 2008 1758 1763 10.1093/PCP/PCN149 18849574
71 Lambour B. Glenz R. Forner C. Krischke M. Mueller M.J. Fekete A. Waller F. Sphingolipid Long-Chain Base Phosphate Degradation Can Be a Rate-Limiting Step in Long-Chain Base Homeostasis Front. Plant Sci. 13 2022 911073 10.3389/FPLS.2022.911073
72 Tsegaye Y. Richardson C.G. Bravo J.E. Mulcahy B.J. Lynch D.V. Markham J.E. Jaworski J.G. Chen M. Cahoon E.B. Dunn T.M. Arabidopsis Mutants Lacking Long Chain Base Phosphate Lyase Are Fumonisin-sensitive and Accumulate Trihydroxy-18:1 Long Chain Base Phosphate J. Biol. Chem. 282 2007 28195 28206 10.1074/JBC.M705074200 17635905
73 Nakagawa N. Kato M. Takahashi Y. Shimazaki K.I. Tamura K. Tokuji Y. Kihara A. Imai H. Degradation of long-chain base 1-phosphate (LCBP) in Arabidopsis: Functional characterization of LCBP phosphatase involved in the dehydration stress response J. Plant Res. 125 2012 439 449 10.1007/S10265-011-0451-9 21910031
74 Gupta P. Roy S. Nandi A.K. MEDEA-interacting protein LONG-CHAIN BASE KINASE 1 promotes pattern-triggered immunity in Arabidopsis thaliana Plant Mol. Biol. 103 2020 173 184 10.1007/S11103-020-00982-4 32100164
75 Systems and Synthetic Biology Laboratory Ssbio/rDMFA (1.0) Zenodo 2024 10.5281/zenodo.12811205
76 DiCiccio T.J. Efron B. Bootstrap confidence intervals Stat. Sci. 11 1996 189 228 10.1214/SS/1032280214
77 Leighty R.W. Antoniewicz M.R. Dynamic metabolic flux analysis (DMFA): A framework for determining fluxes at metabolic non-steady state Metab. Eng. 13 2011 745 755 10.1016/J.YMBEN.2011.09.010 22001431
78 van der Maaten L. Hinton G. Visualizing Data using t-SNE J. Mach. Learn. Res. 9 2008 2579 2605
79 Sander J. Ester M. Kriegel H.P. Xu X. Density-Based Clustering in Spatial Databases: The Algorithm GDBSCAN and Its Applications Data Min. Knowl. Discov. 2 1996 169 194 10.1023/A:1009745219419
80 Isaacs F.J. Hasty J. Cantor C.R. Collins J.J. Prediction and measurement of an autoregulatory genetic module Proc. Natl. Acad. Sci. USA 100 2003 7714 7719 10.1073/pnas.1332628100 12808135
