\chapter{Materials and methods: sphingolipid de novo biosynthesis}
Using a modern control approach, a comparator model is developed to simulate and predict regulatory feedback in a case study metabolic pathway. Specifically, the comparator model is applied to a case study in the effect of single-gene overexpression on sphingolipid \emph{de novo} biosynthesis in a human embryonic kidney (HEK) cell line. From the perspective of interdisciplinary research in particular, the comparator model appeals to the nature of conventional experiment design in the field of biology in terms of traditional case/control studies.
\section{Sphingolipid biology}
Sphingolipids are involved in key eukaryotic cell functions such as membrane structure recognition, signal transduction, intracellular regulation, and cell-cell interaction \cite{merrill-a}. The general structure of a sphingolipid comprises a sphingoid base backbone, e.g., sphinganine, that is modified by addition of long-chain fatty acids, double bonds in the sphingoid base (to make sphingosine), and polar headgroups (Figure 4). When variation in these components is considered, sphingolipids are one of the most complex families of biomolecules \cite{fahy,pruett}.
Long-chain fatty acids and their roles in sphingolipid biosynthesis, metabolism and function are of specific interest, e.g., in ceramide biosynthesis \cite{lahiri}, maintenance of fatty acid levels \cite{choi}, bilayer membrane composition in terms of lipid rafts \cite{megha,garner} and sterol affinity \cite{nyholm}, mitochondrial permeability \cite{novgorodov}, and apoptotic signaling \cite{shabbits}. The metabolism of these sphingolipids is also studied in various cells, e.g., yeast \cite{dickson}, plants \cite{bach}, mouse skin \cite{madison}, and human lung endothelium \cite{medler}.
\emph{De novo} biosynthesis of sphingolipids with varying fatty acid chain lengths, e.g., C16:0 and C26:0, occurs via similar pathways. In particular, they share a common pool of palmitate and sphingoid bases. However, these chain-length variants affect cell functions differently. This suggests that biosynthesis of such chain-length variants is also regulated differently.
Alterations in sphingolipid biosynthesis, storage, and metabolism are implicated in human diseases \cite{kacher}. For example, in synthesis, infantile-onset symptomatic epilepsy is described as a genetic defect in ganglioside GM3 synthase \cite{simpson}; in storage, Tay-Sachs, Sandhoff, and Gaucher disease are characterized by glycosphingolipid accumulation \cite{kolter}; in metabolism, ceramide functions as a tumor suppressor in various cancer cells \cite{ogretmen}. Regulation of sphingolipid levels through \emph{de novo} bio-synthesis is critical to maintain key cell functions, failing which leads to a variety of human diseases. By studying the dynamics and regulation of sphingolipid \emph{de novo} biosynthesis, it may be possible to learn more about onset and development of these sphingolipid diseases.
\clearpage
\begin{figure}[p]%
\begin{center}
\includegraphics[width=3in]{./figures/sph-structure.png}%
\caption[Molecular structure of sphingolipids]
{General molecular structure of sphingolipids with sphingoid base backbone, variable-length fatty acid chains, and headgroups (R).}%
\end{center}
\label{fig:sph-structure}%
\end{figure}
\subsection{Sphingolipid \emph{de novo} biosynthesis}
Sphingolipid \emph{de novo} biosynthesis is a "necessary, but dangerous, pathway" \cite{merrill-a}. \emph{De novo} biosynthesis is necessary because sphingoid bases present in food are mostly degraded in the intestine \cite{vesper}. Regulation of sphingolipid amounts through \emph{de novo} biosynthesis is important because pathway intermediates such as ceramide and sphingomyelin are highly bioactive.
Three of the initial enzymes in the sphingolipid \emph{de novo} biosynthesis pathway (Figure 5) are thought to be particularly important:
\begin{itemize}
\item \emph{serine palmitoyltransferase} (SPT), which catalyzes the initial step of sphingolipid biosynthesis by condensation of L-serine and palmitoyl-coA (Pal-CoA),
\item \emph{(dihydro)ceramide synthase}, which adds the amide-linked fatty acid to form dihydroceramide, and
\item \emph{dihydroceramide desaturase}, which adds the 4,5-trans-double bond of the sphingoid base backbone to form ceramide.
\end{itemize}
\clearpage
\begin{landscape}
\begin{figure}[p]%
\begin{center}
\includegraphics[width=6.5in]{./figures/syn-pathway.png}%
\caption[Sphingolipid \emph{de novo} biosynthesis pathway]
{Sphingolipid \emph{de novo} biosynthesis pathway, with key enzymes in \emph{italics}; uni-/bi-directional arrows indicate irreversible/reversible reactions respectively; pathway intermediates, e.g., ceramide, sphingomyelin, are highly bioactive.}%
\end{center}
\label{fig:syn-pathway}%
\end{figure}
\end{landscape}
\subsubsection{Sphingolipid-omics}
The increasing wealth of quantitative lipidomic data makes the study of sphingolipid biology from a systems perspective promising and challenging at the same time \cite{merrill-c}. The diversity and complexity of sphingolipids require researchers to possess a range of tools to work with these compounds. For example, to study dynamic changes in lipid amounts, computational methods that couple mass spectrometry with statistical algorithms to analyze the vast number of lipid species from cellular extracts are required \cite{forrester}. In addition, the Lipidomics Gateway (www.lipidmaps.org) is a comprehensive website from the LIPID MAPS consortium, the leading group for lipidomics research, which contains standards, data, and tools for researchers interested in lipid biology.
\subsection{Mass spectrometry lipidomics}
Sphingolipids, i.e., sphingoid bases and (dihydro)N-acyl species, are extracted from samples of cultured human embryonic kidney (HEK) cells and quantified by liquid chromatography electrospray ionization tandem mass spectrometry (LC-ESI-MS/MS) as described previously \cite{merrill-b,sullards}. For this case study, measurements from four technical replicates are taken at seven time points in uniform hourly intervals, i.e., at 0, 1, 2, 3, 4, 5, and 6 hours. The community standard is described by the LIPID MAPS consortium, which requires three technical replicates taken at eight time points, i.e., at 0, 0.5, 1, 2, 4, 8, 12, and 24 hours.
The data acquisition protocol for this case study is comparable to the community standard for time-series studies in terms of both (a) number of time points sampled, and (b) number of technical replicates. In addition, compared to the community standard, the data for this case study is also sampled more frequently at the onset of the experiment.
To quantify sphingoid bases labeled with [U-13C]-palmitate, additional multiple reaction monitoring (MRM) pairs corresponding to [M+16] precursor ions and [M+16] product ions are used. To quantify (dihydro)N-acyl sphingolipids labeled with a single [U-13C]-palmitate, MRM pairs corresponding to [M+16] precursor ions and both [M+0] and [M+16] product ions are used, providing discrimination of labeling on the sphingoid base and N-acyl moieties. These [M+16] isotopomers are designated BASE and FA respectively. To quantify (dihydro)N-acyl sphingolipids labeled with two [U-13C]-palmitate molecules, MRM pairs corresponding to [M+32] precursor ions and [M+16] product ions are used; the [M+32] isotopomers are designated DUAL. Unlabeled sphingolipids are designated 12C. The peak areas of unlabeled and labeled sphingolipid isotopomers are integrated, converted to picomoles using the peak area of the internal standard, and normalized to the mg of protein in the extracted sample.
\clearpage
\section{Comparator model}
From engineering control, a \emph{comparator} takes in 2 inputs: one each from the \emph{reference} and the \emph{plant} systems (Figure 6). In general, the reference system exhibits, and establishes, the desired response to an \emph{input}. Under different scenarios, the input could be either known, e.g., a tracking signal, or unknown, e.g., a disturbance to the system. For the purpose of modeling, regardless of whether it is known or unknown, the input is common to both the reference and plant systems. In other words, the reference system is also a \emph{model} for the plant system. Thus, the reference system is fully known while the plant system is unknown (fully or partially) and is therefore the system of interest.
The two inputs to the comparator are the outputs from the reference and plant systems. At each time, the difference between the two given inputs is calculated and used to compute an output by the comparator, which is feedback to the plant system at the next time. This process proceeds iteratively until the plant dynamics converges with the reference dynamics, i.e., when the difference between the outputs from the reference and plant systems is zero (Figure 4).
Thus, the key to convergence between the reference and plant systems lies in how the comparator output is: (a) computed, and (b) implemented, i.e., as actuated feedback to the plant system. As discussed before in Section 1.3.1, the Lyapunov theory of stability provides a direct method to specify how to accomplish the former based on a user-defined stability function.
\clearpage
\begin{landscape}
\begin{figure}[p]%
\begin{center}
\includegraphics[width=6.5in]{./figures/comparator_blkdiag.png}%
\caption[Block diagram: comparator model]
{Block diagram of comparator model: the \emph{comparator} determines the difference in input from the \emph{reference} and \emph{plant} systems, which is feedback to the plant. This process proceeds iteratively such that the plant dynamics converge to the reference dynamics, i.e., the difference between the plant and reference outputs are zero. Both the reference and plant systems received a common \emph{input}.}%
\end{center}
\label{fig:mrac-gen}%
\end{figure}
\end{landscape}
\subsection{Assumptions}
The key assumptions are as follows:
\subsubsection{Treatment effect on system stability}
In this case study, the wild type dynamics is chosen as the \emph{reference} system, which the treated cells, i.e., SPT overexpressed, follow to the same stable steady state. More precisely, pathway dynamics in treated cells is assumed to converge to the wild type over time, i.e., the difference between the two pathways, in terms of measurable metabolite amounts, goes to zero.
Thus, the underlying assumption here is that the effect of the treatment condition, i.e., single-gene (SPT) overexpression, is assumed to be large enough such that pathway dynamics in treated cells can be differentiated from the wild type; however, at the same time, it is small enough such that the treated cells do not converge to a second steady state that is different from the wild type. In other words, homeostasis remains relevant under these experimental conditions. This assumption is supported in theory by the pathway topology \cite{wolf} and in practice from experimental data.
\subsubsection{Determination of pathway differences}
The comparator is a quantitative model of how deviations from the "desired" (stable) state (or response), represented by the wild type in this case, may be handled under homeostasis. The issue is that at any time, while the "desired" state may be observed in wild type cells, only the perturbed state is observed in the treated cells. In other words, without communicating with the wild type in real-time, it is implicitly assumed that the treated cells have knowledge of the "desired" state \emph{a priori}.
Thus, the underlying assumption here is that a "desired" state is already prescribed within the cells, treated or wild type. Such a "desired" state could possibly be genetically predetermined so that if the cell encounters perturbation within certain homeostatic limits, such deviations from the "desired" state could be regulated to restore the system to stability.
\subsection{Relevance to case/control studies in biology}
In this case study, a comparator model is developed and applied to model the effect of single-gene (SPT) overexpression on sphingolipid \emph{de novo} biosynthesis in human embryonic kidney (HEK) cells in terms of pathway regulation under homeostasis. In particular, as a result of homeostasis, treated cells are assumed to regulate sphingolipid \emph{de novo} biosynthesis such that levels of sphingolipid metabolites in these cells approach the wild type over time. Stable isotope 13C-labeled palmitate is added to both treated and wild type cells in excess relative to intracellular sphingolipid metabolites to facilitate tracking the levels of these sphingolipid metabolites by mass spectrometry.
In other words, treated cells are represented as the \emph{plant} system and the wild type as the \emph{reference} system. Sphingolipid precursor, extracellular palmitate is the \emph{input} that is common to both reference and plant systems. Then, by analyzing the dynamics of feedback gain from the \emph{comparator} as a predictive model of homeostasis, interactions between sphingolipid metabolites in the system may be responsible to regulate pathway dynamics (Figure 7).
Thus, the comparator model appeals to, and indeed complements, traditional experiment design in the field of biology in that it provides a ready research paradigm for the analysis of homeostatic mechanisms using conventional case/control experiments, which correspond intuitively to the engineering language of plant/reference systems. Effectively, the proposed comparator model demonstrates that complexities of homeostasis can be reduced by borrowing well-established concepts from control theory to describe equally well-documented phenomena in biology. Consequently, so long as the experimental treatment does not result in responses that are beyond the reference system, the comparator model together with an appropriate theory of stability may be an additional analysis tool that is useful to reverse engineer homeostasis and contribute to research in biochemical robustness.
\clearpage
\begin{landscape}
\begin{figure}[p]%
\begin{center}
\includegraphics[width=6.5in]{./figures/comparator_applied.png}%
\caption[Applied comparator model in a biological case/control study]
{Proposed application of comparator model in a biological case/control study: in this case study on sphingolipid \emph{de novo} biosynthesis, a comparator model from modern engineering control is applied to reverse engineer pathway homeostasis in treated cells (\emph{plant}), i.e., single-gene (SPT) overexpressing cells; the pathway dynamics in the wild type are fully known (\emph{reference}); to both treated and wild type cells, stable isotope 13C-labeled sphingolipid precursor palmitate (\emph{input}) is added in excess relative to intracellular sphingolipid levels, which also faciliate mass spectrometry lipidomics; the key model result is predicted steady-state \emph{feedback} to the treated cells that suggests potential homeostatic pathway interactions as a result of experimental treatment.}%
\end{center}
\label{fig:mrac}%
\end{figure}
\end{landscape}
%\begin{enumerate}
% \item treated cells, i.e., SPT overexpressing cells (plant system), attempt to behave like the wild type (reference system) \-- this assumption is based on the fact that both groups of cells are genetically identical; in other words, aside from SPT overexpression, the cells are clones from the same cell line.
% \item pathway regulation is achieved through adaptive control \-- from a biological perspective, adaptive control \cite{anderson,goodwin} is consistent with current understanding that resources are conserved in metabolism, which is evident from existing studies \cite{friedlander}; from a modeling perspective, compared to other common modes of modern control, there is no need to make further assumptions about cost and performance constraints, for instance, in optimal \cite{bryson,sussman} and robust \cite{van-antwerp-b,van-antwerp-a} control.
%\end{enumerate}
\clearpage
\section{Mathematical representation}\label{algo}
\subsection{State-space model}
There are 10 sphingolipid metabolites of interest in the \emph{de novo} biosynthesis pathway based on the available data. They are:
\begin{enumerate}
\item sphinganine (Sa)
\item sphinganine phosphate (SaP)
\item dihydroceramide (DHCer)
\item dihydroglucosylceramide (DHGC)
\item dihydrosphingomyelin (DHSM)
\item ceramide (Cer)
\item glucosylceramide (GC)
\item sphingomyelin (SM)
\item sphingosine (So)
\item sphingosine phosphate (SoP)
\end{enumerate}
As discussed in Section 2.1.2, extracellular palmitate, a precursor molecule that is labeled using stable isotope 13C, is added to facilitate the measurement of these sphingolipid metabolites using mass spectrometry.
Thus, these sphingolipid metabolites in the \emph{de novo} biosynthesis pathway represent the internal states of the pathway system while the sphingolipid precursor molecule, palmitate, is the input that is common to both wild type and treated cells. Then, the state-space $\vec{x}(t)$ can be defined as the set of intracellular sphingolipid metabolite amounts where
\begin{align}
\vec{x}(t) &= [Sa(t), SaP(t), \nonumber\\
&\qquad DHCer(t), DHGC(t), DHSM(t), \nonumber\\
&\qquad Cer(t), GC(t), SM(t), \nonumber\\
&\qquad So(t), SoP(t)]^{T}
\end{align}
The dynamics of this system can be described in a simplified form using generalized mass action (GMA) \cite{gulberg,lund} and written as a system of linear ordinary differential equations (ODEs) as follows:
\begin{align}
\dot{Sa}(t) &=
k_{in,Sa}PalCoA + k_{SaP,Sa}SaP + k_{DHCer,Sa}DHCer \nonumber\\
&\qquad\qquad - (k_{Sa,SaP} + k_{Sa,DHCer})Sa\\
\dot{SaP}(t)&= k_{Sa,SaP}Sa - k_{SaP,Sa}SaP\\
\dot{DHCer}(t) &=
k_{Sa,DHCer}Sa + k_{DHGC,DHCer}DHGC \nonumber\\
&\qquad + k_{DHSM,DHCer}DHSM + k_{Cer,DHCer}Cer \nonumber\\
&\qquad - (k_{DHCer,Sa} + k_{DHCer,DHGC} + k_{DHCer,DHSM} + k_{DHCer,Cer})DHCer\\
\dot{DHGC}(t)&= k_{DHCer,DHGC}DHCer - k_{DHGC,DHCer}DHGC\\
\dot{DHSM}(t)&= k_{DHCer,DHSM}DHCer - k_{DHSM,DHCer}DHSM\\
\dot{Cer}(t) &=
k_{DHCer,Cer}DHCer + k_{GC,Cer}GC + k_{SM,Cer}SM + k_{So,Cer}So \nonumber\\
&\qquad - (k_{Cer,DHCer} + k_{Cer,GC} + k_{Cer,SM} + k_{Cer,So})Cer\\
\dot{GC}(t)&= k_{Cer,GC}Cer - k_{GC,Cer}GC\\
\dot{SM}(t)&= k_{Cer,SM}Cer - k_{SM,Cer}SM\\
\dot{So}(t)&= k_{Cer,So}Cer + k_{SoP,So}SoP - (k_{So,Cer} + k_{So,SoP})So\\
\dot{SoP}(t)&= k_{So,SoP}So - (k_{SoP,So} + k_{SoP,out})SoP
\end{align}
where $k_{a,b}$ indicates the reaction rate constant for the enzymatic reaction going from $a$ to $b$. In addition, because extracellular palmitate is added in excess relative to intracellular amounts of sphingolipid metabolites, the amount of PalCoA as the system input is assumed to be constant and in excess for the length of the experiment \emph{in vitro} and corresponding simulation \emph{in silico}. On a related note, the Michaelis-Menten model of enzyme reaction kinetics \cite{voet} is not applied here because of the difficulties in estimating the requisite parameters \cite{atkins,currie}.
\subsection{Comparator model}
The comparator model is derived from the Lyapunov theory of stability based on scalar-input to a multi-variable system. As follows in this section, variables and parameters that deal with the \emph{reference} system are denoted by subscript $_{ref}$.
\subsubsection{Plant and reference systems}
The state-space model for the \emph{reference} system, i.e., wild type cells with (assumed) fully known processes, is:
\begin{equation}
\dot{\vec{x}}_{ref}(t) = \mathbf{A}_{ref} \vec{x}_{ref}(t) + \vec{b}_{ref}r(t)
\end{equation}
for the internal states (sphingolipid metabolites) $\vec{x}_{ref}(t)$ of the \emph{reference} system (wild type cells), where dynamics between these states (enzymatic reactions between sphingolipid metabolites) can be described in terms of linear ordinary differential equations (ODEs) with parameters (reaction rate constants $k_{a,b}$) contained in the matrix $\mathbf{A}_{ref}$.
The scalar system \emph{input} (sphingolipid precursor palmitate, which is converted to intracellular palmitoyl-coA PalCoA) $r(t)$ can be scaled and targeted at various internal states by the actuation vector $\vec{b}_{ref}$; in this case, it is involved only in the first enzymatic reaction of the pathway.
The state-space model for the \emph{plant} system, i.e., treated cells with fully or partially unknown processes that are of interest, is:
\begin{equation}
\dot{\vec{x}}(t) = \mathbf{A} \vec{x}(t) + \vec{b}u(t)\\
\end{equation}
for the internal states (sphingolipid metabolites) $\vec{x}(t)$ of the \emph{plant} system (treated cells), where dynamics between these states (enzymatic reactions between sphingolipid metabolites) are described in terms of linear ODEs with parameters (reaction rate constants $k_{a,b}$) contained in the matrix $\mathbf{A}$.
The scalar \emph{feedback} (from PalCoA and other system sphingolipid metabolites) $u(t)$ regulates the dynamics of these various states in the plant through the actuation vector $\vec{b}$. More precisely, the feedback $u(t)$ is the sum of two components such that:
\begin{equation}
u(t) = u_x(t) + u_r(t) = \vec{h}_x(t)^{T}\vec{x}(t) + h_r(t)r(t)
\end{equation}
where $u_x(t)$, $u_r(t)$ are specific to the internal states and system input respectively and are modulated by feedback gain $\vec{h}_x(t)$, $h_r(t)$.
\subsubsection{Derivation from Lyapunov stability}
Let $\vec{e}(t)$ be the difference between \emph{plant} and \emph{reference} systems:
\begin{equation}\label{err}
\vec{e}(t)\equiv \vec{x}(t)-\vec{x}_{ref}(t)
\end{equation}
such that
\begin{align}
\vec{e}(t) &= \vec{x}(t) - \vec{x}_{ref}(t)\nonumber\\
&\Rightarrow\nonumber\\
\dot{\vec{e}}(t) &= \dot{\vec{x}}(t) - \dot{\vec{x}}_{ref}(t)\nonumber\\
&= [\mathbf{A} \vec{x}(t) + \vec{b} u(t)] - [\mathbf{A}_{ref} \vec{x}_{ref}(t) + \vec{b}_{ref} r(t)]\nonumber\\
&= [\mathbf{A} \vec{x}(t) + \vec{b} (\vec{h}_x(t)^{T}\vec{x}(t) + h_r(t)r(t))] - [\mathbf{A}_{ref} \vec{x}_{ref}(t) + \vec{b}_{ref} r(t)]\nonumber\\
&= [\mathbf{A} + \vec{b}\vec{h}_x(t)^{T}]\vec{x}(t)- \mathbf{A}_{ref} \vec{x}_{ref}(t) + [\vec{b}h_r(t) - \vec{b}_{ref}]r(t)\nonumber\\
&= [\mathbf{A} + \vec{b}\vec{h}_{x}^{*T} - \vec{b}\vec{h}_{x}^{*T} + \vec{b}\vec{h}_x(t)^{T}]\vec{x}(t)- \mathbf{A}_{ref} \vec{x}_{ref}(t) + [\vec{b}h_r(t) - \vec{b}h_{r}^{*} + \vec{b}h_{r}^{*} - \vec{b}_{ref}]r(t)\nonumber\\
&= \mathbf{A}_{ref}[\vec{x}(t) - \vec{x}_m(t)] + \vec{b}[\vec{h}_x(t)^{T} - \vec{h}_{x}^{*T}]\vec{x}(t) + \vec{b}[h_r(t) - h_{r}^{*}]r(t)\nonumber\\
&= \mathbf{A}_{ref}\vec{e}(t) + \vec{b}[\Delta\vec{h}_x(t)^{T}\vec{x}(t) + \Delta h_r(t)r(t)]
\end{align}
where
\begin{align}
\mathbf{A} + \vec{b}\vec{h}_{x}^{*T} &= \mathbf{A}_m\nonumber\\
\vec{b}h_{r}^{*} &= \vec{b}_m\nonumber\\
\Delta\vec{h}_x(t) &= \vec{h}_x(t) - \vec{h}_{x}^{*} &\Rightarrow \Delta\dot{\vec{h}}_x(t) = \dot{\vec{h}}_x(t) \nonumber\\
\Delta h_r(t)&= h_r(t) - h_{r}^{*} &\Rightarrow \Delta\dot{h}_r(t) = \dot{h}_r(t)\nonumber
\end{align}
and $\vec{h}_x^{*}$, $h_{r}^{*}$ are the ideal steady-state feedback gains.
Then choose a candidate stability function, i.e., positive definite $V(\cdot)$:
\begin{equation}
V(\vec{e}(t),\Delta\vec{h}_x(t),\Delta h_r(t)) = \vec{e}(t)^{T}\mathbf{P}\vec{e}(t) + \Delta\vec{h}_x(t)^{T}\Gamma_x^{-1}\Delta\vec{h}_x(t) - \gamma_r^{-1}\Delta h_r(t)^{2}
\end{equation}
and take its derivative with respect to time as follows:
\begin{align}
\dot{V}(t) &= \vec{e}(t)^{T}[\mathbf{A}_m^{T}\mathbf{P}+\mathbf{P}\mathbf{A}_m]\vec{e}(t) + 2\vec{e}(t)\mathbf{P}\vec{b}[\Delta\vec{h}_x(t)^{T}\vec{x}(t) + \Delta h_r(t)r(t)] \nonumber\\
&\qquad + 2\Delta\vec{h}_x(t)^{T}\Gamma_x^{-1}\Delta\dot{\vec{h}}_x(t) \nonumber\\
&\qquad + 2\Delta h_r(t)\gamma_r^{-1}\Delta\dot{h}_r(t) \nonumber\\
&= \vec{e}(t)^{T}\mathbf{Q}\vec{e}(t) \nonumber\\
&\qquad + 2\Delta\vec{h}_x(t)^{T}[\vec{x}(t)\vec{e}(t)^{T}\mathbf{P}\vec{b} + \Gamma_x^{-1}\Delta\dot{\vec{h}}_x(t)] \nonumber\\
&\qquad + 2\Delta h_r(t)[r(t)\vec{e}(t)^{T}\mathbf{P}\vec{b} + \gamma_r^{-1}\Delta\dot{h}_r(t)]
\end{align}
Thus, to ensure $\dot{V}(t)$ is negative semi-definite that guarantees asymptotic stability, i.e.,
\begin{equation}
\dot{V}(t) = -\vec{e}(t)^{T}\mathbf{Q}\vec{e}(t) \leq 0
\end{equation}
$\dot{\vec{h}}_x(t)$ and $\dot{h}_r(t)$, which are the rates of change of feedback gain, must be:
\begin{equation}\label{adap-1}
\dot{\vec{h}}_x(t)= -\Gamma_x\vec{x}(t)\vec{e}(t)^{T}\mathbf{P}\vec{b}
\end{equation}
\begin{equation}\label{adap-2}
\dot{h}_r(t)= -\gamma_rr(t)\vec{e}(t)^{T}\mathbf{P}\vec{b}
\end{equation}
where $\Gamma_x$, $\gamma_x$ are also user-defined. Then, $\mathbf{P}$ is positive definite symmetric and satisfies the Lyapunov equation:
\begin{equation}
\mathbf{A}_m^{T}\mathbf{P} + \mathbf{P}\mathbf{A}_m + \mathbf{Q} = \vec{0}
\end{equation}
where $\mathbf{Q}$ is also positive definite symmetric. $\mathbf{P}$ can be solved by specifying $\mathbf{Q}$ in addition to $\mathbf{A}_{ref}$ from equation (1). Finally, the initial conditions $\vec{h}_x(0)$, $h_r(0)$ may be chosen by matching:
\begin{equation}
\mathbf{A} + \vec{b}\vec{h}_{x}(0)^{T} = \mathbf{A}_{ref}
\end{equation}
\begin{equation}
\vec{b}h_{r}(0) = \vec{b}_{ref}
\end{equation}
to approximate the ideal steady-state feedback gains $\vec{h}_x^{*}$, $h_{r}^{*}$.
\clearpage
\section{Implementation}
\subsection{Parameter estimation}
\subsubsection{State-space parameters}
Reaction rate constants $k_{a,b}$ in the state-space model were estimated using a genetic algorithm to optimize a modified least-squares error between \emph{in silico} pathway dynamics and corresponding experimental data \cite{shah}. Only C16:0 DUAL species data from both treated and wild type cell data were used. While additional information may be inferred from each of the labeled/unlabeled combinations of BASE, FA, or 12C species, they are not used here for simplification.
Inputs to the estimation routine include measured amounts of all sphingolipid metabolites determined from 4 technical replicates. Only 5 of 7 data points, from 0 to 4 hours, were used as training data for parameter estimation; the last 2 data points were used to verify estimation results. Furthermore, robustness of \emph{in silico} model results, i.e., dynamics of feedback gain, is also tested by sensitivity analysis of the model results to perturbations in these reaction rate constants, where $k_{perturbed} = k_{original} * {10^{-1}, 10^1}$.
\subsubsection{Comparator parameters}
The system input is extracellular palmitate that is in excess amounts relative to the pathway sphingolipid metabolites. Extracellular palmitate is converted to intracellular palmitoyl-coA (PalCoA). Subsequently, $r(t)$ represents the rate at which PalCoA is converted to sphinganine (Sa), i.e., $k_{in,Sa}PalCoA$, the first enzymatic reaction of the sphingolipid \emph{de novo} biosynthesis pathway. Thus, $r(t)t$ is a boundary condition that is a scalar constant, numerical value 100.
The reference system actuation vector $\vec{b}_{ref}$ supplies only the common system input $r(t)$. The plant system actuation vector $\vec{b}$ supplies the control input in terms of two components: state feedback $u_x(t)$, and system input $u_r(t)$. Numerical values for $\vec{b}$ are specified through iteration, subject to non-negative model results, so as to minimize the least-squares error between \emph{in silico} model results and \emph{in vitro} experimental data.
$\Gamma_x$, $\gamma_r$, and $\mathbf{Q}$ are parameters of the candidate stability function, subject to conditions of the Lyapunov theory of stability. These parameters affect the rate at which the plant system approaches (stable) steady state, and in practice, i.e., for controller design, may be tuned to achieve specific performance objectives. Here, as a first attempt to model homeostasis in a highly regulated metabolic pathway, these parameters are selected non-specifically.
$\Gamma_x$ and $\gamma_r$ correspond to the internal states of the pathway system $x(t)$ and the system input $r(t)$ respectively. Starting from an initial value of the identity matrix $\mathbf{I}$ for $\Gamma_x$ and unity for $\gamma_r$, these parameters are tuned iteratively, by one order of magnitude to $0.1$, to alleviate the stiffness problem for numerical integration \cite{quo-a}.
$\mathbf{Q}$ is specified as the identity matrix $\mathbf{I}$. Then, based on the Lyapunov equation, $\mathbf{A}_{ref}^{T}\mathbf{P} + \mathbf{P}\mathbf{A}_{ref} + \mathbf{Q} = \vec{0}$, $\mathbf{P}$ can be solved in MATLAB using the built-in function \texttt{lyap()} once $\mathbf{Q}$ and $\mathbf{A}_{ref}$ are also specified.
$h_r(0)$ and $\vec{h}_x(0)$ represent initial conditions for feedback gains that correspond to the system input $r(t)$ and state-space $x(t)$. The numerical values are specified as $1$ for $h_r(0)$ and $\vec{0}$ for $\vec{h}_x(0)$. Based on the candidate stability function subject to the Lyapunov theory of stability, the rate of change of these feedback gains is a function of the state-space difference $\vec{e}(t)$ between the reference and plant systems. In theory, initial conditions for the feedback gains affect only the \emph{rate} of convergence; however in practice, if these initial conditions are poorly chosen, then convergence may not be achieved at all.
\subsection{Modeling error}
The modeling error is measured in terms of the root-sum-square (RSS) quantity between \emph{in silico} model results and \emph{in vitro} experimental data, which is defined as:
\begin{equation}\label{rss}
RSS = \sqrt{\sum^{t_f}_{t=t_0}(x_{in silico,t}-x_{in vitro,t})^2}
\end{equation}
where $x_{in silico,t}$ represents \emph{in silico} model results at time $t$ and $x_{in vitro,t}$ represents \emph{in vitro} experimental data at time $t$, where $t = {1, 2, 3, 4, 5, 6}$.
\subsection{Numerical details}
The aforementioned mathematical representation is implemented in MATLAB using the built-in function \texttt{ode15s()} as follows:
\begin{itemize}
\item $\mathbf{A},\mathbf{A}_{ref},\Gamma_x,\mathbf{P},\mathbf{Q}\in\Re^{10x10}$
\item $\vec{x}(t),\vec{x}_{ref}(t),\vec{e}(t),\vec{h}_x(t),\vec{b}, \vec{b}_{ref}\in\Re^{10}$
\item $u(t), r(t),h_r(t),\gamma_r\in\Re$
\end{itemize}
The proposed comparator model is first applied to data on the C16:0 sphingolipid \emph{de novo} biosynthesis pathway. Then, to verify the generality of the proposed model, it is subsequently applied to data on the C26:0 sphingolipid pathway. While the two sphingolipid pathways are similar in terms of \emph{de novo} biosynthesis, they are independent systems because:
\begin{itemize}
\item different enzymes are responsible to catalyze enzymatic reactions in each pathway system -- a \emph{family} of (dihydro)ceramide synthases (CerS1, CerS2, …, CerS6) are responsible for dihydroceramide (DHCer) synthesis from sphinganine (Sa); in particular, CerS3 and CerS4 produce DHCer with very long chain fatty acids that are greater than or equal to 20 carbon atoms long,
\item C16:0 and C26:0 sphingolipid metabolites are quantified using mass spectrometry independently,
\item parameters of the state-space models, i.e., reaction rate constants, are estimated independently.
\end{itemize}
In reality, the underlying biological picture is likely to be much more complex. For instance, there may be cross-signaling events between the two pathways simply because sphingolipid metabolites are localized in the same regions on cellular membranes. Nonetheless in this case study, because such cross-signaling events in \emph{de novo} biosynthesis between sphingolipid metabolites that are significantly different in fatty acid chain length is not yet well documented, modeling such events is not yet an issue.
\subsubsection{C16:0 pathway}
Parameters of the C16:0 pathway state-space model, i.e., reaction rate constants $k_{a,b}$, are available in Appendix A. All other parameters are defined as follows:
\begin{align*}
\texttt{input: } & r(t) = 100 ( = k_{in,Sa}PalCoA) \\
\texttt{actuation: } & \vec{b}_{ref} = [1, 0, \ldots, 0]^{T},\\
& \vec{b} = [10, -0.1, -10, 10, 10, -0.1, -0.1, -0.1, -0.1, 10]^{T}\\
\texttt{adaptation: } & \mathbf{Q} = \mathbf{I}, \qquad \Gamma_x = 0.1\mathbf{I}, \qquad \gamma_r = 0.1\\
\texttt{initial conditions: } & h_r(0) = 1,\\
& \vec{h}_x(0)=[0, \ldots, 0]^{T}\\
\end{align*}
\subsubsection{C26:0 pathway}
Parameters of the C26:0 pathway state-space model, i.e., reaction rate constants $k_{a,b}$, are available in Appendix B. All other parameters are defined as follows:
\begin{align*}
\texttt{input: } & r(t) = 100 ( = k_{in,Sa}PalCoA) \\
\texttt{actuation: } & \vec{b}_{ref} = [1, 0, \ldots, 0]^{T},\\
& \vec{b} = [10, -0.1, -10, 10, 10, -0.1, -0.1, -0.1, -0.1, 10]^{T}\\
%& \vec{b} = [-1, 0, 1, 1, -15, 0, 1, 0, -1, 0; 1, 0, \ldots, 0]^{T}\\
\texttt{adaptation: } & \mathbf{Q} = \mathbf{I}, \qquad \Gamma_x = 0.1\mathbf{I}, \qquad \gamma_r = 0.1\\
\texttt{initial conditions: } & h_r(0) = 1,\\
& \vec{h}_x(0)=[0, \ldots, 0]^{T}\\
\end{align*}