\chapter{Results: pathway dynamics and feedback}
In this chapter, results of the proposed comparator model to reverse engineer homeostasis in a highly regulated metabolic pathway are reported as follows:
\begin{itemize}
\item to test (a) the efficacy of the proposed comparator model in capturing observed pathway dynamics, and (b) the robustness of predicted steady-state pathway feedback in response to simulated errors in state-space parameter estimation,
\item to interpret and verify predicted steady-state pathway feedback in the context of homeostasis in C16:0 sphingolipid \emph{de novo} biosynthesis, by comparison to literature as well as independent data on measured changes in sphingolipid metabolite levels in response to 4HPR dosage in both treated and wild type cells, and
\item to demonstrate the generality of the proposed modeling approach, by applying the comparator model to a different metabolic system using independent data on C26:0 sphingolipid metabolites.
\end{itemize}
Here, to facilitate discussion of the model results, the 10 sphingolipid metabolites of interest in the \emph{de novo} biosynthesis pathway can be grouped as follows:
\begin{itemize}
\item Sa* \-- Sa and SaP,
\item So* \-- So and SoP,
\item DH* \-- DHCer, DHGC, and DHSM,
\item Cer* \-- Cer, GC, and SM.
\end{itemize}
where Sa*, So* are sphingoid bases, DH* are dihydrosphingolipids, and Cer* are sphingolipids.
\clearpage
\section{Model efficacy and robustness}
\subsection{Pathway dynamics}
In this subsection, the model results are reported in terms of: (a) state-space parameter fold change between treated and wild type cells, (b) root-sum-square (RSS) error between \emph{in silico} model results and \emph{in vitro} experimental data in wild type cells, i.e., \emph{reference} system, (b) RSS error between \emph{in silico} model results and \emph{in vitro} experimental data in treated cells, i.e., \emph{plant} system, \emph{without} comparator, and (d) RSS error between \emph{in silico} model results and \emph{in vitro} experimental data in treated cells \emph{with} comparator.
\subsubsection{State-space parameter fold change between treated and wild type cells}
Each state-space parameter represents the reaction rate constant of a specific enzymatic reaction in the pathway. Appendix A lists the estimated state-space parameters for the C16:0 sphingolipid biosynthesis pathway. As such, parameter fold change between comparable systems is often used as a preliminary indicator of enzymatic reactions that may proceed at noticeably different rates in one system than the other, i.e., in treated cells compared to the wild type, which may be biologically interesting.
Fold change in reaction rate constants between treated and wild type cells is reported in Table 2. Reaction rate constants with observed $log_{10}$ fold change greater than magnitude "$1$" may be considered to be "biologically significant". In addition, based on these values, Figure 10 shows the $log_{10}$ fold change of these reaction rate constants between treated and wild type cells. Fold change is "significant" in 9 of the 18 reaction rate constants in all.
The 9 reaction rate constants with potentially "biologically significant" fold change between treated and wild type cells are:
\begin{description}
\item[fold change $<10^{-2}$] $k_{Sa,DHCer}$, $k_{GC,Cer}$,
\item[$10^{-2}<$ fold change $<10^{-1}$] $k_{DHCer,Sa}$, $k_{DHCer,DHGC}$, $k_{DHCer,Cer}$, $k_{Cer,So}$, $k_{So,Cer}$, and
\item[$10^1<$ fold change $<10^2$] $k_{Sa,SaP}$, $k_{DHGC,DHCer}$.
\end{description}
\clearpage
\begin{table}[p]
\caption{Parameter fold change: C16:0 reaction rate constants}
\begin{center}
\begin{tabular}{lr@{.}l}
\hline
\hline
% $k_{in,Sa}$ & 98&9476 \\
\noalign{\smallskip}
$k_{SaP,Sa}$ & -0&1811\\
\noalign{\smallskip}
$k_{DHCer,Sa}$ & -1&2927\\
\noalign{\smallskip}
$k_{Sa,SaP}$ & 1&3986\\
\noalign{\smallskip}
$k_{Sa,DHCer}$ & -2&8536\\
\noalign{\smallskip}
$k_{DHGC,DHCer}$ & 1&3260\\
\noalign{\smallskip}
$k_{DHSM,DHCer}$ & -0&1277\\
\noalign{\smallskip}
\hline
\hline
\noalign{\smallskip}
$k_{Cer,DHCer}$ & 0&5080\\
\noalign{\smallskip}
$k_{DHCer,DHGC}$ & -1&5468\\
\noalign{\smallskip}
$k_{DHCer,DHSM}$ & 0&4930\\
\noalign{\smallskip}
$k_{DHCer,Cer}$ & -1&1099\\
\noalign{\smallskip}
$k_{GC,Cer}$ & -3&7469\\
\noalign{\smallskip}
$k_{SM,Cer}$ & 0&1851\\
\noalign{\smallskip}
\hline
\hline
\noalign{\smallskip}
$k_{So,Cer}$ & -1&6749\\
\noalign{\smallskip}
$k_{Cer,GC}$ & 0&0072\\
\noalign{\smallskip}
$k_{Cer,SM}$ & -0&6839\\
\noalign{\smallskip}
$k_{Cer,So}$ & -1&2724\\
\noalign{\smallskip}
$k_{SoP,So}$ & -0&1414\\
\noalign{\smallskip}
$k_{So,SoP}$ & 1&4213\\
\noalign{\smallskip}
\hline
\hline
\end{tabular}
\label{table:p-fold-change}
\end{center}
\end{table}
\clearpage
\begin{landscape}
\begin{figure}[p]%
\begin{center}
\includegraphics[width=6.5in]{./figures/k16-comp_v02.png}%
\caption[Parameter fold change: C16:0 reaction rate constants]
{Parameter fold change between treated cells (\emph{plant}) over wild type (\emph{reference}): colors indicate range of fold change in terms of logarithmic ratio $log_{10}$ of estimated reaction rate constants of C16:0 sphingolipid \emph{de novo} biosynthesis pathway.}%
\end{center}
\label{fig:k16foldchange}%
\end{figure}
\end{landscape}
\subsubsection{Wild type cells (\emph{reference} system)}
In the following series of time-series graphs of model results, the scatterplots represent median value and range of 4 technical replicates (where available) of \emph{in vitro} experimental data in discrete time ($t = {0, 1, 2, ..., 6}$) and broken lines represent \emph{in silico} model results in terms of sphingolipid metabolite amounts.
Figure 8 shows the experimental data and model results for pathway dynamics in wild type cells in terms of sphingolipid metabolite amounts using the generalized mass action (GMA) state-space model. Experimental data was not available for Sa*.
For So*, the median data suggests that So increases quickly at first until time $t = 2h$ and approaches a steady state around $t = 4h$. At the same time, measurements for So are less consistent compared to other molecules and have a much larger range. The median values for SoP over time shows that SoP remains low throughout. This observation is supported by its range of measurements, which is consistently small except for at time $t = 6h$.
For DH*, technical replicates were not available and the measurements fluctuate significantly throughout the course of experiment. As such, it is not possible to estimate apparent trends in how the amounts of these molecules change over time.
For Cer*, experimental data is highly consistent across all replicates as seen from the small range of measurements. Thus, the data shows that Cer increases quickly at first till $t = 1h$ and approaches a steady state around $t = 4h$; GC increases very slightly throughout $t = 0h$ to $t = 6h$, and SM increases more quickly before $t = 3h$ compared to after $t = 4h$.
Comparing the model results with experimental data qualitatively, the state-space model of pathway dynamics based on GMA is a good fit to the experimental data for So*, Cer, and GC (4 metabolite species), where they fall clearly and consistently within the range of measurements for these metabolites. For SM, the model result reflects similar trends, i.e., an increase over time but at a much larger rate \emph{in silico} than is measured \emph{in vitro}. However for DH*, the model results are a poor fit to the experimental data, although even to the trained eye, it may not be straightforward to discern an apparent trend in pathway dynamics for these sphingolipid metabolites given the limited data.
\clearpage
\begin{figure}[p]%
\begin{center}
\includegraphics[width=4.5in]{./figures/c16-0_wt.png}%
\caption[Wild type C16:0 sphingolipid metabolites]
{Experimental data and model results for C16:0 sphingolipid metabolites in wild type cells (\emph{reference} system): abscissa \-- time (hours), ordinate -- experimentally normalized sphingolipid metabolite amounts (pmol/mg protein); scatterplots represent \emph{in vitro} experimental data (median and range indicated where available), broken lines represent \emph{in silico} model results.}%
\end{center}
\label{fig:c16-0_wt}%
\end{figure}
\clearpage
\subsubsection{Treated cells (\emph{plant} system) without comparator}
Figure 9 shows the experimental data and model results of pathway dynamics in treated cells in terms of sphingolipid metabolite amounts using the generalized mass action (GMA) state-space model \emph{without} the comparator.
For Sa*, the median data suggests that Sa increases very quickly at first until time $t = 2h$, then decreases almost immediately towards a steady state at $t = 3h$. However, the rate of the initial increase may be overestimated given that the measurements for Sa at $t = 1, 2h$ are spread over a larger range compared to measurements at subsequent times. The data for SaP is more consistent throughout the course of the experiment, which shows that SaP increases only slightly over time.
For So*, the data suggests that So increases continually throughout the course of experiment, although the range of measurements for So is much larger at times $t = 5, 6h$ compared to the preceding time. The data for SoP is similar to that for SaP and is consistent throughout, which shows that SoP also increases only slightly over time.
Technical replicates were not available for DH*. However, the available data does suggest gradual changes over time. DHCer appears to increase quickly from $t = 1h$ to $t = 2h$, then decreases more gradually from $t = 3h$ to $t = 5h$. DHGC appears to remain unchanged throughout the course of experiment. DHSM appears to increase quickly from $t = 1h$ to $t = 2h$ before approaching a steady state subsequently.
For Cer*, the experimental data is highly consistent across all replicates for GC and SM, as can be seen from the small range of measurements, but not for Cer where the range of measurements for Cer is much larger at times $t = 2, 3h$ compared to at other times. The data suggests that Cer increases quickly from $t = 1h$ to $t = 2h$, then decreases gradually from $t = 2h$ to $t = 5h$. DHGC increases only very slightly till $t = 6h$. SM increases uniformly and at a slightly higher rate over the same period of time.
Comparing the model results with experimental data qualitatively, \emph{in silico} sphingolipid metabolite amounts based on the GMA state-space model \emph{without} comparator are a good fit to the data for 5 metabolite species (SaP, SoP, DHGC, GC, and SM), where they fall clearly and consistently within the range of measurements for these metabolites. However, for 4 metabolite species (Sa, So, DHCer, DHSM, and Cer), the model results consistently underestimate the experimental data. Furthermore, for 2 species (DHCer and Cer), the model results do not reflect changes in the experimental data sufficiently, in particular where rapid increases in the measured data from $t = 2h$ to $t = 3h$ followed by more gradual decreases from $t = 1h$ to $t = 5h$ are not captured in the state-space model.
\clearpage
\begin{figure}[p]%
\begin{center}
\includegraphics[width=4.5in]{./figures/c16-0_spt.png}%
\caption[Treated cells without comparator: C16:0 sphingolipid metabolites]
{Experimental data and model results for C16:0 sphingolipid metabolites in treated cells (\emph{plant} system) \emph{without} comparator: abscissa \-- time (hours), ordinate -- experimentally normalized sphingolipid metabolite amounts (pmol/mg protein); scatterplots represent \emph{in vitro} experimental data (median and range indicated where available), broken lines represent \emph{in silico} model results.}%
\end{center}
\label{fig:c16-0_spt}%
\end{figure}
\clearpage
\subsubsection{Root-sum-square (RSS) error in wild type and treated cells state-space model}
Compared to model results for the wild type (Figure 8), model results for the treated cells (Figure 9) are not as good a fit to the experimental data. In particular, the model results for the treated cells are observed to have the poorest during intervals where significant changes are recorded in the measurement data. In terms of root-sum-square (RSS) error, the goodness-of-fit of the GMA state-space model of both wild type and treated cells \emph{without} comparator is quantified and reported in Table 3. The RSS error is also equivalent to the L2-norm between model results and experimental data taken at discrete times, i.e., where $t = {0, 1, 2, ..., 6}$. Taken together, these observations show that although the same algorithm was used for parameter estimation, i.e., a modified genetic algorithm described in Section 2.4.1, model results based on these estimated state-space parameters still differ in accuracy of matching experimental data for different datasets.
\clearpage
\begin{table}[p]
\caption{State-space error: C16:0 wild type and treated cells \emph{without} comparator}
\begin{center}
\begin{tabular}{@{}lcc@{}}
\hline
\noalign{\smallskip}
& Wild type & Treated cells $^{\natural}$\\
\noalign{\smallskip}
\hline
\hline
\noalign{\smallskip}
C16:0 sphingolipids &&\\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
Sa & 0.0002 & 61.4536\\
\noalign{\smallskip}
\noalign{\smallskip}
SaP & 0.0004 & 5.5469\\
\noalign{\smallskip}
\noalign{\smallskip}
DHCer & 6.9245 & 327.3739\\
\noalign{\smallskip}
\noalign{\smallskip}
DHGC & 16.8494 & 3.0832\\
\noalign{\smallskip}
\noalign{\smallskip}
DHSM & 6.0871 & 171.2229\\
\noalign{\smallskip}
\noalign{\smallskip}
Cer & 4.7858 & 213.9073\\
\noalign{\smallskip}
\noalign{\smallskip}
GC & 9.1085 & 41.2140\\
\noalign{\smallskip}
\noalign{\smallskip}
SM & 88.0422 & 19.6755\\
\noalign{\smallskip}
\noalign{\smallskip}
So & 20.2075 & 99.8462\\
\noalign{\smallskip}
\noalign{\smallskip}
SoP & 2.4847 & 9.4587\\
\noalign{\smallskip}
\hline
\end{tabular}
\label{table:wt-v-trt}
\end{center}
\begin{flushleft}
Accuracy of state-space GMA model in terms of root-sum-square (RSS) error of \emph{in silico} model results compared to \emph{in vitro} experimental data in wild type and $^{\natural}$ treated cells \emph{without} comparator.
%$^{\flat}$ in terms of root-sum-square (RSS) error, i.e., square of the sum of differences between numerical results (computed sphingolipid amounts) and median data (from experiment) across discrete time: 0, 1, 2, \ldots, 6h\\
%$^{\sharp}$\emph{percentage change} over RSS difference in treated cells without MRAC in parentheses.
\end{flushleft}
\end{table}
\clearpage
\subsubsection{Treated cells (\emph{plant} system) with comparator}
Figure 10 shows the experimental data and model results of pathway dynamics in treated cells in terms of sphingolipid metabolite amounts using the generalized mass action (GMA) state-space model \emph{with} the comparator. In particular, the model results are a good fit to the experimental data for 4 metabolite species (So, SoP, DHGC, and DHSM). The model results underestimate the measurement data for 1 metabolite species (Sa), and overestimate the data for 3 metabolite species (SaP, GC, and SM). For 2 metabolite species (DHCer and Cer), the model results do not capture the transient trend of the sphingolipid pathway dynamics from time $t = 1h$ to $t = 5h$.
In addition, to illustrate the effect of the comparator model to simulate state-space pathway dynamics in the treated cells, Table 5 shows the quantitative accuracy of the model results to experimental data in terms of comparison between the treated cells \emph{with} and \emph{without} comparator. It is noteworthy that the quantitative comparison shows that model accuracy in simulating state-space pathway dynamics is comparable in both the existing approach using GMA \emph{without} comparator as in the proposed approach \emph{with} comparator.
\clearpage
\begin{figure}[p]%
\begin{center}
\includegraphics[width=5in]{./figures_v02/c16-0_spt_mrac.png}%
\caption[Treated cells with comparator: C16:0 sphingolipid metabolites]
{Experimental data and model results for C16:0 sphingolipid metabolites in treated cells (\emph{plant} system) \emph{with} comparator: abscissa \-- time (hours), ordinate -- experimentally normalized sphingolipid metabolite amounts (pmol/mg protein); scatterplots represent \emph{in vitro} experimental data (median and range indicated where available), broken lines represent \emph{in silico} model results.}%
\end{center}
\label{fig:c16-0_spt_mrac}%
\end{figure}
\clearpage
\begin{table}[p]
\caption{State-space error: C16:0 treated cells \emph{with} and \emph{without} comparator}
\begin{center}
\begin{tabular}{@{}lcc@{}}
\hline
\noalign{\smallskip}
& \emph{Without} comparator & \emph{With} comparator\\
\noalign{\smallskip}
\hline
\hline
\noalign{\smallskip}
Treated cells only &&\\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
Sa & 61.4536 & 70.1442\\
\noalign{\smallskip}
\noalign{\smallskip}
SaP & 5.7469 & 8.2915\\
\noalign{\smallskip}
\noalign{\smallskip}
DHCer & 327.3739 & 346.0480\\
\noalign{\smallskip}
\noalign{\smallskip}
DHGC & 3.0832 & 2.9421\\
\noalign{\smallskip}
\noalign{\smallskip}
DHSM & 171.2229 & 176.7805\\
\noalign{\smallskip}
\noalign{\smallskip}
Cer & 213.9073 & 227.8527\\
\noalign{\smallskip}
\noalign{\smallskip}
GC & 41.2140 & 46.6592\\
\noalign{\smallskip}
\noalign{\smallskip}
SM & 19.6755 & 11.2155\\
\noalign{\smallskip}
\noalign{\smallskip}
So & 99.8462 & 29.9782\\
\noalign{\smallskip}
\noalign{\smallskip}
SoP & 9.4587 & 4.8419\\
\noalign{\smallskip}
\hline
\end{tabular}
\label{table:unctrl-v-ctrl}
\end{center}
\begin{flushleft}
%$^{\flat}$ in terms of root-sum-square (RSS) error, i.e., square of the sum of differences between numerical results (computed sphingolipid amounts) and median data (from experiment) across discrete time: 0, 1, 2, \ldots, 6h\\
%$^{\sharp}$\emph{percentage change} over RSS difference in treated cells without MRAC in parentheses.
\end{flushleft}
\end{table}
%\clearpage
%\subsubsection{Model efficacy}
%In terms of model efficacy, Figure 12 shows the tracking error of both controlled and uncontrolled plants to the reference system in terms of root-sum-square error (RSS) between simulation results for the treated (plant) and wild type (reference) cells. Of 10 variables, 9 (except SM) have decreased tracking error in the controlled plant compared to the uncontrolled plant. In other words, the tracking error of most variables in the controlled plant is improved over the uncontrolled plant.
%
%\begin{landscape}
%\begin{figure}[p]%
% \begin{center}
% \includegraphics[width=5.5in]{./figures_v02/rss-efficacy_org.png}%
% \caption[Model efficacy: tracking error in controlled/uncontrolled plant]
% {Tracking error in terms of root-sum-square (RSS) for metabolites in treated cells: uncontrolled plant, i.e., mass action only (blue) and controlled plant, i.e., with adaptive feedback (red).}%
% \end{center}
% \label{fig:C16:0-rss_error}%
%\end{figure}
%\end{landscape}
%\clearpage
%\subsection{Robustness of feedback}
%Because simulation results are directly affected by parameter noise, the robustness of the controlled plant to improve tracking error is further assessed by perturbing the kinetic parameters systematically. Specifically, each parameter is perturbed separately over a range of 2 orders of magnitude, i.e., $k_{perturbed} = k_{original} * {10^{-1}, 10^1}$.
%
%Figure 13 shows the mean and standard deviation of the tracking error in response to these perturbations. There are two key observations that can be made from these results: $(a)$ the tracking error is smaller in the controlled plant for 9 of 10 state variables (pathway metabolites) compared to the uncontrolled plant, and $(b)$ this relationship is robust under perturbation, i.e., when the kinetic parameters are modified systematically.
%
%\clearpage
%\begin{landscape}
%\begin{figure}[p]%
% \begin{center}
% \includegraphics[width=5.5in]{./figures_v02/rss-efficacy_prtb.png}%
% \caption[Model efficacy: sensitivity analysis of tracking error]
% {Mean tracking error in terms of root-sum-square (RSS) as a result of parameter sensitivity analysis /-- kinetic parameters in treated cells were perturbed systematically over 2 orders of magnitude, i.e., $k_{perturbed} = k_{original} * {10^{-1}, 10^1}$: uncontrolled plant, i.e., mass action only (blue), and controlled plant, i.e., with adaptive feedback (red); error bars indicate standard deviation; $p < 0.05$ for Sa, SaP, and DHGC using two-tailed $t$-test}%
% \end{center}
% \label{fig:C16:0-eff-prtb}%
%\end{figure}
%\end{landscape}
\clearpage
\subsection{Steady-state feedback}
Figures 12 and 13 shows the key model outcome, i.e., predicted steady-state pathway feedback in terms of input and aggregate state feedback as well as individual sphingolipid metabolite feedback based on the proposed comparator model of homeostasis in C16:0 sphingolipid \emph{de novo} biosynthesis in response to single-gene (SPT) overexpression in HEK cells.
Figure 12 shows the predicted steady-state pathway feedback in treated cells in terms of 2 components: (a) system input and sphingolipid precursor, palmitoyl-CoA (PalCoA), $u_r(t)$, and (b) aggregate state feedback from pathway sphingolipid metabolites $u_x(t)$. At steady-state, the magnitude of the feedback suggests the contribution of each component to homeostasis, i.e., a larger magnitude indicate greater contribution, and sign (+/-) indicates either enhancement (positive feedback) or inhibition (negative feedback) to the pathway. The range indicated at discrete time, i.e., $t = {1, 2, ..., 6}$, is the result of sensitivity analysis of these model results in response to potential errors in state-space parameter estimation (as discussed in Section 2.4.1).
From the onset, i.e., at time $t = 0h$, $u_r(t)$ oscillates for a very brief period of time but with large peak-to-peak amplitude, then quickly decreases from the zero initial condition to settle at steady-state (approximately $-20$). On the other hand, throughout the course of simulation, $u_x(t)$ increases steadily from the zero initial condition to settle at steady-state (approximately $+30$). At steady-state, the magnitude of feedback in each of these components is similar but in different directions, i.e., the different signs suggest that input feedback $u_r(t)$ is inhibited while aggregate state feedback $u_x(t)$ is enhanced as a result of homeostasis.
\clearpage
\begin{figure}[p]%
\begin{center}
\includegraphics[width=4.5in]{./figures_v02/c16-0_spt_ctrl_v02.png}%
\caption[Predicted steady-state feedback: C16:0 pathway input and aggregate states]
{Predicted steady-state feedback for C16:0 sphingolipid metabolites pathway in treated cells in terms of: (a) system input and sphingolipid precursor, palmitoyl-coA (PalCoA), $u_r(t)$ (blue), and (b) aggregate state feedback from pathway sphingolipid metabolites $u_x(t)$ (red) - x-axis: time (hours), y-axis: magnitude (dimensionless), error bars: standard deviation from state-space parameter sensitivity analysis}%
\end{center}
\label{fig:c16-0_spt_ctrl}%
\end{figure}
\clearpage
Figure 13 shows the predicted steady-state feedback gain from individual pathway sphingolipid metabolites $h_x(t)$. These feedback gain are dimensionless; however, at steady-state, the magnitude and sign of these feedback values may be interpreted to suggest differential activity in pathway homeostasis. More precisely, the larger the magnitude, the greater the feedback; furthermore, the sign (+/-) of the gain values indicate either enhancement (positive feedback) or inhibition (negative feedback) to the pathway. In addition, the range indicated at discrete time, i.e., $t = {1, 2, ..., 6}$, is the result of sensitivity analysis of these model results in response to potential errors in state-space parameter estimation (as discussed in Section 2.4.1).
Sa, SaP, So, and SoP are sphingoid bases; DHCer, DHGC, DHSM, Cer, GC, and SM are sphingolipids. Thus, from Figure 13:
\begin{description}
\item[\emph{(top left)}] Feedback gain for sphinganine (Sa) decreases quickly from the zero initial condition to settle at steady-state (approximately $-1$); there is no noticeable change in feedback gain for sphinganine-phosphate (SaP).
\item[\emph{(top right)}] There is no noticeable change in feedback gain for sphingosine (So); however, feedback gain decreases quickly from the zero initial condition to settle at steady-state (approximately $-1$) for sphingosine-phosphate (SoP).
\item[\emph{(bottom left)}] Feedback gain for dyhydroceramide (DHCer) increases quickly from the zero initial condition to settle at steady-state (approximately $+1$); on the other hand, it decreases quickly from the zero initial condition to settle at steady-state (approximately $-1$) for dihydroglucosylceramide (DHGC) and dihydrosphingomyelin (DHSM).
\item[\emph{(bottom right)}] Feedback gain for ceramide (Cer) decreases gradually from the zero initial condition and approaches $-0.07$ at time $t = 6h$; there is no noticeable change from the zero initial condition for glucosylceramide (GC) and sphingomyelin (SM).
\end{description}
\clearpage
\begin{figure}[p]%
\begin{center}
\includegraphics[width=4.5in]{./figures_v02/c16-0_spt_fdbk.png}%
\caption[Predicted steady-state feedback: individual C16:0 sphingolipid metabolites]
{Predicted steady-state feedback for individual C16:0 sphingolipid metabolites in treated cells - x-axis: time (hours), y-axis: magnitude (dimensionless), error bars: standard deviation from state-space parameter sensitivity analysis}%
\end{center}
\label{fig:c16-0_spt_fdbk}%
\end{figure}
\clearpage
\section{Biological interpretation and verification}
\subsection{Pathway input and aggregate state feedback}
From Figure 12, these model results suggest that in response to single-gene (SPT) overexpression in treated cells, and when sphingolipid precursor PalCoA is added to the cells in excess amounts relative to intracellular sphingolipid metabolite levels, the metabolic effect of the pathway input (PalCoA) at steady-state is to inhibit sphingolipid \emph{de novo} biosynthesis. On the other hand, the collective metabolic effect of pathway sphingolipid metabolites at steady-state is to enhance \emph{de novo} biosynthesis. This observation suggests that, in this case of serine pamitolytransferase (SPT) upregulation (as a result of single-gene overexpression) where SPT is the enzyme that catalyzes for the initial pathway reaction of condensation between PalCoA and L-serine to form sphinganine, the excess input and reaction substrate PalCoA also inhibits the pathway dynamics.
This implication is intriguing in that in terms of controller design, such a regulatory strategy is also very specifically tuned to the type of system input. From a biochemical perspective, given that there may be only a limited number of admissible pathway inputs, such a regulatory strategy may be a viable option to maintain pathway homeostasis. For instance, there could be molecular mechanisms to slow the rate of extracellular palmitate uptake \cite{glatz,luiken,mcarthur}, reduce the rate of conversion from palmitate to PalCoA, or even increase the rate of other enzymatic reactions that require PalCoA as a substrate. Subsequently, given that PalCoA is an integral component of sphingolipid metabolites at least in terms of molecular structure, the implication that regulation of the sphingolipid \emph{de novo} biosynthesis pathway may be highly dependent on intracellular PalCoA levels should not be too surprising.
The possibility that levels of intracellular PalCoA may be moderated as a result of SPT overexpression could also be interpreted as a result of experimental procedure, i.e., where only data on dual-chain isotope-labeled sphingolipid metabolites was used to develop this model (see Section 2.4.1 on parameter estimation). It could be that in reality, regardless of isotope-labeling, sphingolipid \emph{de novo} biosynthesis and turnover is increased in SPT overexpressed cells, which leads to increased palmitate that is subsequently converted back to PalCoA so that it is available to re-enter the pathway. Subsequently, because only data on dual-labeled sphingolipids is used for parameter estimation, it may appear that the amount of 13C-labeled PalCoA is moderated relative to the levels of all intracellular PalCoA, labeled or otherwise. This scenario may be tested experimentally by improving the tracking of various pools of labeled and unlabeled palmitate and PalCoA using mass spectrometry.
\clearpage
\subsection{Individual sphingolipid metabolite feedback}
Steady-state feedback gain in individual sphingolipid metabolites may be interpreted as follow: relative to other sphingolipid metabolites such as So and Cer, steady-state feedback gain from Sa, SoP, DHCer, DHGC, and DHSM is significantly different from the zero initial condition. This suggests that these sphingolipid metabolites may contribute more actively in pathway homeostasis. Furthermore, these particular sphingolipid metabolites also contribute to cell signaling. In particular, Sa and SoP are also cell signaling molecules that are reported to activate growth stimulation: Sa reverses growth inhibition \cite{solomon} while SoP activates growth stimulation \cite{spiegel-a,spiegel-b}. In addition, So and Cer are also reported to activate growth inhibition \cite{cuvillier,merrill-d,woodcock}.
Taken together in the context of pathway homeostasis, the similarities and differences in steady-state feedback gain for these sphingolipid signaling metabolites, which activate opposing responses, suggests that the experimental treatment, i.e., SPT overexpression at these levels, leads to growth stimulation in HEK cells \cite{wei}. Clinically, SPT overexpression is also implicated in cancer cell metastasis in human tumors \cite{carton}.
Thus, in general, model results and interpretation from Figures 12 and 13 are consistent with current knowledge of sphingolipid biology (summarized in Table 5). The robustness of these model results, i.e., predicted steady-state feedback, in response to potential errors in parameter estimation underscores the utility and promise the proposed comparator model to simulate and predict homeostatic mechanisms in a highly regulated metabolic pathway. Furthermore, these model results and present interpretations raise rather new and perhaps interesting questions with regard to more specific aspects of homeostasis in sphingolipid \emph{de novo} biosynthesis pathway that were not previously studied.
\clearpage
\begin{landscape}
\begin{table}[p]
\caption{Sphingolipid metabolite signaling: reported effects on cell growth}
\begin{center}
\begin{tabular}{lp{5cm}p{5.8cm}}
\hline
\noalign{\smallskip}
\noalign{\smallskip}
Sphingolipid metabolite & Reported effects on cell growth & Predicted steady-state feedback\\
\noalign{\smallskip}
\noalign{\smallskip}
\hline
\hline
\noalign{\smallskip}
\noalign{\smallskip}
Sphinganine (Sa) & Growth stimulation \cite{solomon} & Strong inhibition\\
\noalign{\smallskip}
\noalign{\smallskip}
Sphinganine phosphate (SaP) & & No change\\
\noalign{\smallskip}
\noalign{\smallskip}
Sphingosine (So) & Growth inhibition \cite{cuvillier,merrill-d,woodcock} & Weak inhibition\\
\noalign{\smallskip}
\noalign{\smallskip}
Sphingosine phosphate (SoP) & Growth stimulation \cite{spiegel-a,spiegel-b} & Strong inhibition\\
\noalign{\smallskip}
\noalign{\smallskip}
Dihydroceramide (DHCer) & Not available & Strong enhancement\\
\noalign{\smallskip}
\noalign{\smallskip}
Dihydroglucosylceramide (DHGC) & Not available & Strong inhibition\\
\noalign{\smallskip}
\noalign{\smallskip}
Dihydrosphingomyelin (DHSM) & Not available & Strong inhibition\\
\noalign{\smallskip}
\noalign{\smallskip}
Ceramide (Cer) & Growth inhibition \cite{cuvillier,merrill-d,woodcock} & Weak inhibition\\
\noalign{\smallskip}
\noalign{\smallskip}
Glucosylceramide (GC) & Not available & Weak enhancement\\
\noalign{\smallskip}
\noalign{\smallskip}
Sphingomyelin (SM) & Not available & Weak enhancement\\
\noalign{\smallskip}
\noalign{\smallskip}
\hline
\end{tabular}
\label{table:fdbk-interpretation}
\end{center}
\begin{flushleft}
\end{flushleft}
\end{table}
\end{landscape}
\clearpage
\subsection{Additional testing with $1 {\mu}M$ 4HPR}
4HPR (\emph{4-hydroxyphenylretinamide}) is a synthetic retinoid that is under clinical evaluation as a therapeutic agent in a variety of cancers. It is reported to restrain tumor growth by inducing apoptosis. However, the required effective dosage to achieve toxicity differs significantly in SPT overexpressing cells and wild type HEK cells. In particular, a higher dose of 4HPR is required in SPT overexpressing cells compared to wild type HEK cells to achieve observable, or effective, cell toxicity (Figure 14).
In particular, although the mode of action, i.e., molecular mechanisms, of 4HPR has been linked to ceramide metabolism (Figure 15), its systemic effects are not fully understood. With respect to the three key enzymes involved in de novo sphingolipid synthesis, 4HPR is reported to enhance the activity of (a) SPT, i.e., serine palmitoyl-transferase that catalyzes the reaction of serine and palmitoyl-coA to form sphinganine (Sa), (b) (DH)CerS, i.e., (dihydro)ceramide synthase that catalyzes the formation of (dihydro)ceramide (DHCer) from sphinganine (Sa), and inhibits the activity of (c) DES, i.e., desaturase that catalyzes the formation of ceramide (Cer) from (dihydro)ceramide (DHCer). Subsequently, from a modeling perspective, the effect of 4HPR on sphingolipid \emph{de novo} biosynthesis can be interpreted as a change in state-space parameters because the reaction rate constants are altered. In particular, 4HPR treatment is similar to SPT overexpression to the extent that SPT is overexpressed in both treatments.
\clearpage
\begin{figure}[p]%
\begin{center}
\includegraphics[width=4.5in]{./figures/4HPR_dosage.png}%
\caption[Effect of 4HPR on SPT and HEK cell viability]
{Effect of 4HPR on SPT and HEK cell viability -- larger dose of 4HPR is required to achieve effective toxicity in SPT cells compared to wild type HEK cells: y-axis --number of viable cells (\% of cells without 4HPR, experimental data), x-axis -- 4HPR dosage}%
\end{center}
\end{figure}
\clearpage
\begin{figure}[p]%
\begin{center}
\includegraphics[width=4.5in]{./figures/4HPR_sphingometabolism.png}%
\caption[Action of 4HPR on sphingolipid metabolism]
{Action of 4HPR on sphingolipid metabolism: 4HPR enhances the activity of enzymes SPT, (DH)CerS, and inhibits DES}%
\end{center}
\end{figure}
\clearpage
As a result of 4HPR treatment, from the known mechanisms of 4HPR (Figure 15), intracellular levels of Sa and DHCer are expected to increase while Cer levels are expected to decrease in HEK cells. Experimental data in Figure 16 illustrates the measured changes in the ratio of intracellular sphingolipid metabolite levels DH*:Cer* at time $t = 24h$ after 4HPR treatment, which is as expected in HEK cells based on the current understanding of 4HPR mechanisms on sphingolipid metabolism (Figure 15). However, the same expected changes in sphingolipid metabolite levels in SPT cells are not observed.
Because 4HPR treatment is similar to SPT overexpression to the extent that SPT activity is enhanced in both treatments, model results in terms of predicted steady-state feedback is verified in part by comparing Figures 13 and 16, where the predicted steady-state feedback (Figure 13) is quantitatively consistent with the observed changes in ratios of sphingolipid metabolite levels at time $t = 24h$ based on the available experimental data (Figure 16). In other words, quantitative agreement of the predicted steady-state feedback with additional independent experimental data suggests that the model predictions can be used to identify, without prior knowledge, which pathway metabolites may be responsible for differences in SPT vs. HEK cell response to 4HPR treatment, i.e., which metabolites may be controllable or are more involved in pathway homeostasis.
From a second perspective (that is not necessarily exclusive from the first), the predicted steady-state feedback specific to modeling homeostasis as a result of SPT overexpression could also contribute to a more comprehensive study to understand the metabolic effects of 4HPR by deconstructing its known mechanisms independently. In this case, modeling the effect of SPT overexpression, and its predicted steady-state feedback from available data, is the first step in testing and modeling individual mechanisms of 4HPR separately before these effects are combined. Consequently, quantitative agreement between the model results, i.e., predicted steady-state feedback, with independent 4HPR response data should lend confidence that the proposed approach using a comparator model is viable, given that the proposed approach and comparator model is built only on experimental data from SPT overexpression that represents only one of three reported mechanisms of 4HPR. Subsequently, future work may apply the proposed approach similarly but on other data from (DH)CerS overexpression and DES underexpression separately in order to complete the 4HPR response model.
Thus, the proposed approach using a comparator model to reverse engineer homeostasis in metabolic pathways is useful in providing (a) a numerical approach to predict quantitative differences in sphingolipid metabolite levels observed between SPT and HEK cells as a result of specific experimental treatments, and (b) an analytical paradigm, in terms of homeostasis and control, to try and answer \emph{why} these differences are observed.
\clearpage
\begin{figure}[p]%
\begin{center}
\includegraphics[width=4.5in]{./figures/4HPR_sphingolevels.png}%
\caption[Effect of $1 {\mu}M$ 4HPR on DHCer/Cer levels (experimental data)]
{Effect of $1 {\mu}M$ 4HPR on DHCer/Cer levels (experimental data): DHCer:Cer ratios are increased in HEK cells, which is consistent with expectation of increased DHCer and decreased Cer levels; changes in SPT ratios are consistent with predicted steady-state feedback (Figure 13)}%
\end{center}
\end{figure}
%C26:0 pathway
\clearpage
\section{Generality of proposed approach}
\subsection{Application to C26:0 sphingolipid metabolites}
To demonstrate the generality of the proposed approach, the comparator model is applied to a second, independent dataset of C26:0 sphingolipid metabolites. Figures 17 and 18 show the experimental data and model results, in terms of sphingolipid metabolite amounts, for wild type and treated cells based on a similar GMA state-space model of pathway dynamics. For the wild type (Figure 17), no data were available for Sa*, and technical replicates were not available for DH*. For treated cells (Figure 18), technical replicates were not available for DH*.
C16:0 and C26:0 sphingnolipid metabolites have the same pathway topology for \emph{de novo} biosynthesis, which also draws on a common pool of sphingolipid precursor PalCoA as the initial pathway reaction substrate. However, experimental data shows that the pathway dynamics of C26:0 sphingolipid metabolites is very different from the C16:0 sphingolipid metabolites. In both wild type and treated cells, similar amounts of sphingoid bases, So*, are observed for both the C16:0 and C26:0 sphingolipid pathways, but slightly less dihydrosphingolipids, DH*, and significantly less sphingolipids, Cer*, are detected in the C26:0 pathway compared to the C16:0 pathway. Because homeostasis is driven primarily by changes in system dynamics, differences in pathway dynamics between both systems from experimental data should be captured in the model results as well.
In Figure 17 for wild type cells, for So*, experimental data show that the amount of So increases gradually from $t = 0h$ and settles to steady-state at $t = 4h$. The amount of SoP remains low over the course of experiment. For DH*, noticeable changes are observed in DHCer over the course of experiment but without any obvious discernible trend, while amounts of DHGC and DHSM remain low throughout. For Cer*, the experimental data show that no dual-labeled Cer, GC, or SM were detected. The model results do not fit well for most metabolites, where they tend to overestimate the experimental data.
In Figure 18 for treated cells, for Sa*, the experimental data shows that Sa increases quickly from $t = 0h$ to $t = 2h$, then decreases to settle at steady-state by $t = 6h$. SaP increases gradually over the course of experiment. For So*, both So and SoP increase steadily throughout the course of experiment, but So increases at a much faster rate. For DH*, the experimental data show that no dual-labeled DHCer, DHGC, or DHSM were detected. For Cer*, Cer and SM increase slightly after $t = 2h$ and $t = 5h$ correspondingly while no GC is detected. Qualitatively, the model results are a good fit with the experimental data for 5 metabolite species (SaP, SoP, DHGC, DHSM, and GC), fair for 4 metabolite species (Sa, So, Cer, and SM), and poor for 1 metabolite species (DHCer).
Figure 19 shows the experimental data and model results for C26:0 sphingolipid metabolites in treated cells with the comparator model. Qualitatively, the model results are a good fit for 5 metabolite species (SaP, SoP, DHGC, DHSM, and GC), fair for 3 metabolite species (Sa, So, and Cer), and poor for 2 metabolite species (DHCer and SM). Table 6 shows a quantitative measure of the accuracy of these model results in matching experimental data at discrete time in terms of RSS error. With regard to the GMA state-space model of pathway dynamics (top 2 rows), the model results are a better fit to the experimental data for the wild type than treated cells. Comparing the model results for treated cells with and without the comparator model (bottom 2 rows), goodness-of-fit differs most significantly for DH*, with an increase in RSS error of $19.4\%$ with the comparator model, while RSS error is decreased for the other groups of sphingolipid metabolites. Overall, RSS error increases by $18.9\%$ with the comparator model, where the increase in RSS error for DH* is most likely the primary factor for the difference.
\clearpage
\begin{figure}[p]%
\begin{center}
\includegraphics[width=4.5in]{./figures/c26-0_wt.png}%
\caption[Wild type C26:0 sphingolipid metabolites]
{Experimental data and model results for C26:0 sphingolipid metabolites in wild type cells (\emph{reference} system): abscissa \-- time (hours), ordinate -- experimentally normalized sphingolipid metabolite amounts (pmol/mg protein); scatterplots represent \emph{in vitro} experimental data (median and range indicated where available), broken lines represent \emph{in silico} model results.}%
\end{center}
\label{fig:c26-0_wt}%
\end{figure}
\clearpage
\begin{figure}[p]%
\begin{center}
\includegraphics[width=4.5in]{./figures/c26-0_spt.png}%
\caption[Treated cells without comparator: C26:0 sphingolipid metabolites]
{Experimental data and model results for C26:0 sphingolipid metabolites in treated cells (\emph{plant} system) \emph{without} comparator: abscissa \-- time (hours), ordinate -- experimentally normalized sphingolipid metabolite amounts (pmol/mg protein); scatterplots represent \emph{in vitro} experimental data (median and range indicated where available), broken lines represent \emph{in silico} model results.}%
\end{center}
\label{fig:c26-0_spt}%
\end{figure}
\clearpage
\begin{figure}[p]%
\begin{center}
\includegraphics[width=4.5in]{./figures/c26-0_spt_mrac.png}%
\caption[Treated cells with comparator: C26:0 sphingolipid metabolites]
{Experimental data and model results for C26:0 sphingolipid metabolites in treated cells (\emph{plant} system) \emph{with} comparator: abscissa \-- time (hours), ordinate -- experimentally normalized sphingolipid metabolite amounts (pmol/mg protein); scatterplots represent \emph{in vitro} experimental data (median and range indicated where available), broken lines represent \emph{in silico} model results.}%
\end{center}
\label{fig:c26-0_spt_mrac}%
\end{figure}
\clearpage
\begin{table}[p]
\caption{State-space error: C26:0 sphingolipid metabolites}
\begin{center}
\begin{tabular}{@{}llccccc@{}}
\hline
\noalign{\smallskip}
& & Sa* & DH* & Cer* & So* & all\\
\noalign{\smallskip}
\hline
\hline
\noalign{\smallskip}
C26:0 & & & & & &\\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
Wild type & & n/a & 13.78 & 0.29 & 458.23 & 458.43\\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
Treated cells & \emph{without} comparator & 54.10 & 860.87 & 0.86 & 104.56 & 868.88\\
\noalign{\smallskip}
$^{\sharp}$ & \emph{with} comparator & 48.88 & 1027.90 & 0.82 & 95.07 & 1033.40\\
&& \emph{(-9.6)} & \emph{(+19.4)} & \emph{(-4.7)} & \emph{(-9.1)} & \emph{(+18.9)}\\
\noalign{\smallskip}
\hline
\end{tabular}
\label{table6}
\end{center}
\begin{flushleft}
%$^{\flat}$based on root-sum-square (RSS), calculated using experimental data (median values) across all discrete time, i.e., 0, 1, 2, \ldots, 6h;\\
$^{\sharp}$ percentage change over RSS error in treated cells \emph{without} comparator
\end{flushleft}
\end{table}
\clearpage
\subsubsection{Steady-state feedback}
Figures 20 and 21 shows the key model outcome, i.e., predicted steady-state pathway feedback in terms of input and aggregate state feedback as well as individual sphingolipid metabolite feedback based on the proposed comparator model of homeostasis in C26:0 sphingolipid \emph{de novo} biosynthesis in response to single-gene (SPT) overexpression in HEK cells. Because the C26:0 sphingolipid data and biology is different from the C16:0 metabolites, it is expected that the model results, i.e., predicted steady-state feedback, are different from the latter case.
Figure 20 shows the predicted steady-state pathway feedback in treated cells in terms of 2 components: (a) system input and sphingolipid precursor, palmitoyl-CoA (PalCoA), $u_r(t)$, and (b) aggregate state feedback from pathway sphingolipid metabolites $u_x(t)$. At steady-state, the magnitude of the feedback suggests the contribution of each component to homeostasis, i.e., a larger magnitude indicate greater contribution, and sign (+/-) indicates either enhancement (positive feedback) or inhibition (negative feedback) to the pathway.
Both components exhibit a brief period of transient oscillation from $t = 0h$ to $t = 1h$ with similar frequency. However, the peak-to-peak amplitude of oscillation is significantly larger for $u_r(t)$ at $\sim270$ compared to $u_x(t)$ at $\sim140$. In addition, $u_r(t)$ settles to steady-state at $\sim+130$ while $u_x(t)$ settles to steady-state at $\sim0$. Compared to C16:0 sphingolipids, the dynamics of feedback control in C26:0 sphingolipids are clearly different in terms of period and amplitude of transient oscillations. At the same time, these feedback dynamics are consistent with designing a control policy that utilizes the least resources for pathway regulation, which is similar to the case for C16:0 sphingolipids.
Figure 21 shows the predicted steady-state feedback gain from individual pathway sphingolipid metabolites $h_x(t)$. These feedback gain are dimensionless; however, at steady-state, the magnitude and sign of these feedback values may be interpreted to suggest differential activity in pathway homeostasis. More precisely, the larger the magnitude, the greater the feedback; furthermore, the sign (+/-) of the gain values indicate either enhancement (positive feedback) or inhibition (negative feedback) to the pathway.
The steady-state feedback gain is non-zero for 5 metabolite species (PalCoA, So, DHGC, DHSM, and GC), but is unchanged from the zero initial condition for the other sphingolipid metabolites. Thus, Figures 19 and 20 show that the predicted steady-state feedback gain for C26:0 sphingolipid metabolites are different from the C16:0 sphingolipid metabolites, which suggests a different set of homeostatic interactions may be responsible for pathway regulation in this system. On the other hand, predicted steady-state feedback due to system input palmitate $u_r(t)$ for both C16:0 and C26:0 sphingolipid metabolites are similar in that PalCoA is suggested to have a self-regulatory role, where excess input PalCoA appears also to inhibit the dynamics of sphingolipid \emph{de novo} biosynthesis.
From a modeling perspective, the significance of these model results is that the proposed approach, using a comparator model to reverse engineer homeostasis in a highly regulated metabolic pathway, can be applied to other datasets in general. Because of a lack of prior knowledge, the biological interpretation and verification of these model results are not discussed in detail here. Furthermore, the model results are consistent with the understanding that homeostasis is driven primarily by pathway dynamics, as opposed to pathway topology, which can be derived directly from time-series experimental data.
\clearpage
\begin{figure}[p]%
\begin{center}
\includegraphics[width=4.5in]{./figures/c26-0_spt_mrac_sysreg.png}%
\caption[Predicted steady-state feedback: C26:0 pathway input and aggregate states]
{Predicted steady-state feedback for C26:0 sphingolipid metabolites pathway in treated cells in terms of: (a) system input and sphingolipid precursor, palmitoyl-coA (PalCoA), $u_r(t)$ (blue), and (b) aggregate state feedback from pathway sphingolipid metabolites $u_x(t)$ (red) - x-axis: time (hours), y-axis: magnitude (dimensionless)}%
\end{center}
\label{fig:c26-0_spt_mrac_sysreg}%
\end{figure}
\clearpage
\begin{figure}[p]%
\begin{center}
\includegraphics[width=4.5in]{./figures/c26-0_spt_mrac_idvfdb.png}%
\caption[Predicted steady-state feedback: individual C26:0 sphingolipid metabolites]
{Predicted steady-state feedback for individual C26:0 sphingolipid metabolites in treated cells - x-axis: time (hours), y-axis: magnitude (dimensionless)}%
\end{center}
\label{fig:c26-0_spt_mrac_idvfdb}%
\end{figure}
%\clearpage
%\section{Comparing C16:0 and C26:0 sphingolipids}
%
%Subsequently, the proposed MRAC approach suggests markedly different regulation dynamics for the C26:0 pathway compared to C16:0. Relatively large, quick bursts in the systematic control response are observed for C26:0 while the equivalent C16:0 response is sustained over longer periods of time. Similar events are also observed for feedback from individual molecules in C26:0, indicating that steady state is reached earlier for C26:0 than for C16:0.
%
\clearpage
\section{Summary}
In summary, \emph{in silico} results of the proposed approach using a comparator model to reverse engineer homeostasis in a highly regulated metabolic pathway have shown:
\begin{itemize}
\item (a) the efficacy of the proposed comparator model in capturing observed pathway dynamics, and (b) robustness of predicted steady-state pathway feedback in response to simulated errors in state-space parameter estimation,
\item the validity of predicted steady-state pathway feedback in the context of homeostasis in C16:0 sphingolipid \emph{de novo} biosynthesis, by comparison to literature as well as independent data on measured changes in sphingolipid metabolite levels in response to 4HPR dosage in both treated and wild type cells, and
\item the generality of the proposed modeling approach, by applying the comparator model to a different metabolic system using independent data on C26:0 sphingolipid metabolites.
\end{itemize}
Using C16:0 sphingolipid \emph{de novo} biosynthesis as a case study, the proposed approach to reverse engineer pathway homeostasis using a comparator model resulted in the prediction of steady-state pathway feedback, which is qualitatively and quantitatively verified with literature as well as additional independent experimental data. Furthermore, the comparator model is also applied to a separate independent dataset on C26:0 sphingolipids, where differences in model results are also consistent with differences in experimental data between the 2 datasets. As an initial study in "the use of existing techniques in well-developed areas of control theory to analyze problems of interest to biologists" \cite{sontag}, the outcomes of this case study on C16:0 sphingolipid metabolites have also been published in \cite{quo,quo-b}.