%\documentclass[10pt]{article}
\documentclass[10pt]{report}

\usepackage{amsmath,amsfonts}
\usepackage{url} 
\usepackage{fullpage} 
\usepackage[numbers,sort&compress,sectionbib]{natbib} % \cite,\citet,\citep,\ldots
\usepackage{graphicx,psfrag}
\usepackage{ifthen}        % \ifthenelse,\boolean,\newboolean,\setboolean
\usepackage{subfigure}     % \subfigure[a]{\includegraphics\ldots\label{fi:\ldots}}
%\usepackage{sidecaption}   % \sidecaption (to be placed inside figure env.

\usepackage[dvips,colorlinks=true,linkcolor=blue,breaklinks=true]{hyperref}
\usepackage[norefs,nocites]{refcheck}  % options:norefs,nocites,msgs,chkunlbld
%\usepackage[nocites]{refcheck}  % options:norefs,nocites,msgs,chkunlbld
\usepackage{mathptmx}      % slightly more compressed font. 
\usepackage{verbatim}
\usepackage{makeidx}       % to make a keyword index
%\usepackage{showidx}      % prints index entries in the left margin (debug)
\usepackage{chemarr}

%%%%%%%%%%%%%% source code formating for latex %%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{listings}
%%%% Matlab language
\lstset{language=Matlab}
\lstset{basicstyle=\tt}
\lstset{commentstyle={\footnotesize\sl}}
\lstset{frame=none}
%\lstset{frame=shadowbox}
%\lstset{columns={[l]fixed}}  % character alignment
\lstset{columns={[l]fixed},basewidth={.55em,.4em}}  % character alignment
%\lstset{columns={[l]fullflexible},basewidth={.6em,.45em}}  % character alignment
\lstdefinelanguage[extended]{Matlab}[]{Matlab}{morekeywords={linprog,quadprog,ones,optimset,ode23s,expand}}

%%%% StochDynTools dialect of Matlab, based on Matlab without dialect
\lstdefinelanguage[StochDynTools]{Matlab}[extended]{Matlab}{morekeywords={readNet,closureDynamics,mylatex,subsParameters,momentDynamics,momentClosure,momentEquilibrium,plotCMoments,getDistribution,net2stochKit,centerMoments}}
%%%% language for .net files
\lstdefinelanguage{net}{basicstyle=\tt,morecomment=[l]{\%},keywords={species,parameters,reactions,substitutions,stochastic,constant,deterministic,rate,derivatives}}

%%%% set default language
\lstset{language=[StochDynTools]Matlab}
%\lstset{language=Matlab}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\usepackage[draft,fancythm,fancybb,morse]{jphmacros2e} 
\usepackage[fancythm,fancybb,morse]{jphmacros2e} 

%% Macros & options for this Document

%\newcommand{\SDT}{\texttt{StochDynTools}}

\newcommand{\codesize}{\footnotesize}

\newcommand{\toidx}[1]{\index{\lstinline{#1}}}%\addcontentsline{toc}{subsubsection}{\lstinline{#1}}}

\newenvironment{command}[1]{\toidx{#1}\addcontentsline{toc}{subsection}{\lstinline{#1}}\paragraph*{\lstinline[basicstyle=\large,columns={[l]flexible}]{#1}}~\\\noindent\rule{\textwidth}{2pt}\\\vspace{-3ex}\codesize}{\vspace{-3ex}\rule{\textwidth}{1pt}\medskip\noindent}

\makeindex

\allowdisplaybreaks

%% Start of Document

\title{\sc \SDT \\[1em]\Large A \matlab{} Toolbox to Compute Moment
  Dynamics\\ for Stochastic Networks of Bio-Bhemical Reactions}

\author{\jph}

\date{July 1, 2009}

\begin{document}                        \maketitle

\begin{abstract}
  This document provides a guide to use the \SDT{} \matlab{} toolbox
  \cite{HespanhaDec06}.  The core of this library is a set of
  functions that compute the approximate moment dynamics for a network
  of chemical reactions.  Examples are provided to illustrate the use
  of the library.

\vfill

\textsc{Attention:} This toolbox is still in a very early stage of
development. It has been posted online mostly for the use of people
that I collaborate with. However, anyone is welcome to try it.
\emph{Please let me know if you have problems or questions about it.}
Your help in improving this package will be greatly appreciated!

\end{abstract}


\tableofcontents

\chapter{Quick start}\label{se:quick}

The \SDT{} toolbox provides functions to compute the approximate moment
dynamics for a network of chemical reactions. The approximate moment
dynamics are based on moment closure techniques that approximate
higher-order moments as static functions of lower-order moments
\cite{HespanhaSinghSep05,SinghHespanhaDec06,UribeVerghese06}. This
section shows by example the use of the \SDT{} toolbox.

\medskip

Networks of chemical reactions are specified by a \texttt{.net} file.
The syntax of a \texttt{.net} file is somewhat self-explanatory.  For
example, consider the Michaelis-Menten mechanism for enzyme kinetics
\begin{align*}
  S + EF &\xrightleftharpoons[D_C]{K_C} ES \xrightarrow{K_P} EF + P
\end{align*}
where $EF$ denotes the (free) enzyme, $S$ the substrate, $ES$ the
enzyme-substrate complex, and $P$ the product. This network of three
elementary chemical reactions can be described by the following
\texttt{.net} file:
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5.25in}
\begin{lstlisting}[language=net,columns={[l]flexible},basewidth={.6em,.45em}]
species:
        EF stochastic  2;   % free enzyme 
        ES stochastic  1;   % enzyme-substrate complex
        S  stochastic  9;   % free substrate
        P  stochastic  1;   % product
parameters:
        K_C "k_C" =  1;    % rate constant for  S + EF -> ES
        D_C "d_C" = 20;    % rate constant for  ES -> S + EF
        K_P "k_P" =.05;    % rate constant for  ES -> P + EF
reactions:
        % reversible reaction
        rate = K_C*S*EF;   {S,EF,ES} > {S-1,EF-1,ES+1};   % S + EF -> ES 
        rate = D_C*ES;     {S,EF,ES} > {S+1,EF+1,ES-1};   % ES -> S + EF
        % product creation reaction
        rate = K_P*ES;     {P,EF,ES} > {P+1,EF+1,ES-1};   % ES -> P + EF
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}  
In this network, the reversible reaction $S + EF \xrightleftharpoons{}
ES$ typically occurs much more often than the production of $P$. This
means that the total number of substrate molecules $St= S + ES$ and
the total number of enzyme molecules $Et=EF+ES$ have much slower
dynamics that those of the free number of substrate molecules $S$ and
the free number of enzymes $EF$. In fact, for the above network $Et$ is
constant and therefore its stochasticity can be ignored. Because of
this, it is more convenient to keep track of the ``slow'' variables
$Et$ and $St$, instead of the ``fast'' variables $ES$ and $S$. Note
that the number of enzyme-substrate complex molecules can be obtained
from $ES=Et-EF$ and the number of free substrate molecules can be
obtained from $S=St-ES=St-Et+EF$.  This approach was followed in the
following \texttt{.net} file, in which $EF$ is the only ``fast''
variable.
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5.25in}
\begin{lstlisting}[language=net,columns={[l]flexible},basewidth={.6em,.45em}]
% Filename: enzyme_with_Et_St.che
species:
        EF "e_f" stochastic   2;  % free enzyme                     (fast variable)
        Et "e_t" constant     3;  % total enzymes = free+compound   (constant)
        St "s_t" stochastic  10; % total substrate = free+compound (slow variable)
        P  "p"    stochastic   1;  % product                         (slow variable)
parameters:
        K_C "k_C" =  1;        % rate constant for  S + EF -> ES
        D_C "d_C" = 20;        % rate constant for  ES -> S + EF
        K_P "k_P" =.05;        % rate constant for  ES -> P + EF
reactions:
        % fast reversible reaction
        rate = K_C*(St-Et+EF)*EF;  {EF} > {EF-1};                 % S + EF -> ES
        rate = D_C*(Et-EF);          {EF} > {EF+1};               % ES -> S + EF
        % slow product creation reaction
        rate = K_P*(Et-EF);          {St,P,EF} > {St-1,P+1,EF+1}; % ES -> P + EF
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}  
The \lstinline{readNet()} command of \SDT{} can be used to read the
\lstinline{enzyme_with_Et_St.net} file above:
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5.25in}
\begin{lstlisting}
net=readNet('enzyme_with_Et_St.net');
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}  
The following \SDT{} command then computes the first-order approximate
moment dynamics, which is roughly the \emph{deterministic chemical rate}
equation\footnote{The main difference being that reactions of the form
  $A+A\xrightarrow{} B$ gives rise to terms of the form $\pm k A(A-1)$,
  instead of $\pm k A^2$.}:
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5.25in}
\begin{lstlisting}
mdyn1=closureDynamics(net,1,'derMatch');
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}  
One could produce \emph{\LaTeX{} code} for these approximate moment dynamics
using the following commands
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5.25in}
\begin{lstlisting}
dotMu=expand(mdyn1.approxDotMu.sym);
latexmacros=mylatex()
latexMu=mylatex(mdyn1.Mu.sym,mdyn1.texrules)
latexDotMu=mylatex(dotMu,mdyn1.texrules)
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}  
where, if one were to uncomment the second line, the parameters would
be replaced by their numerical values. The \LaTeX{} code in the
variables \lstinline{latexMu} and \lstinline{latexDotMu} requires the macros
in the variable \lstinline{latexmacros}. Copying the content of
\lstinline{latexMu} and \lstinline{latexDotMu} into a \LaTeX{} equation
environment one obtains
\begin{align*}
\Deriv{}{t}
\matt{\E[e_f]\\
 \E[s_t]\\
 \E[p]}
=
\matt{-\E[e_f]d_C-k_P \E[e_f]+\E[e_f]k_C e_t+d_C e_t+k_P e_t-k_C{\E[e_f]}^{2}-k_C \E[e_f]\E[s_t]\\
 k_P \E[e_f]-k_P e_t\\
-k_P \E[e_f]+k_P e_t}
\end{align*}
Once can compute the \emph{equilibrium point} of the system above
using the following \matlab{} commands
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5.25in}
\begin{lstlisting}
muEq=momentEquilibrium(net,mdyn1);
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}  
The following \SDT{} command then computes the \emph{second-order
  approximate moment dynamics} using derivative matching
(cf.~Appendix~\ref{se:derivative-matching}) and produces an m-file
\lstinline{fun.m} that computes the derivatives of the uncentered moments
of order up to 2:
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5.25in}
\begin{lstlisting}
mdyn2=closureDynamics(net,2,'derMatch','fun');
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}  
The following \matlab{} commands solve the approximate moment dynamics
and plot the evolution of the \emph{uncentered moments}:
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5.25in}
\begin{lstlisting}
Tmax=300;                      % maximum simulation time
[t,mu]=ode23s(@(t,x)fun(x),[0,Tmax],mdyn2.Mu.x0);
plot(t,mu)
legend(fun(),'location','best','interpreter','none')
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}  
One could instead plot the evolution of the \emph{centered moments} of
\lstinline{E}, \lstinline{St}, and \lstinline{Et} using the following \SDT{} command
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5.25in}
\begin{lstlisting}
h=plotCMoments(net,mdyn2,t,mu,'EF','b-','St','g-','Et','k-');
legend(h(1:3:end),'E[EF]\pm{}Std[EF]','E[St]\pm{}Std[St]','E[Et]','location','best');
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}  
One could also plot the \emph{distribution} of \lstinline{St} using the
following \SDT{} command
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5.25in}
\begin{lstlisting}
x=0:.01:4;
pdf=getDistribution(net,mdyn2,mu(end,:),'St',x,'lognormal');
plot(x,pdf)
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}  
The outputs of the previous plot commands are shown in
Figure~\ref{fig:plots}.
\begin{figure}[htbp]
  \centering
  \hfill
  \subfigure[uncentered moments]{\includegraphics[width=.3\textwidth]{uncentered.eps}}
  \hfill
  \subfigure[centered moments]{\includegraphics[width=.3\textwidth]{centered.eps}}
  \subfigure[distribution of St]{\includegraphics[width=.3\textwidth]{pdf.eps}}
  \caption{Output of the plot commands in the example in
    Section~\ref{se:quick}.}\label{fig:plots}
\end{figure}

One could also run Monte Carlo simulations of the system using the
following \SDT{} command
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5.25in}
\begin{lstlisting}
[Q,b,c,s,x0]=quadPropensities(net);
Q=double(subsParameters(net,Q));      % replace parameters
b=double(subsParameters(net,b));      % by their numerical values
c=double(subsParameters(net,c));
s=double(subsParameters(net,s));
nMC=1000;                             % number of Monte Carlo runs
Ts=(0:1:300)';                        % simulation times of interest
[X,Xmean,Xstd]=sampledSSA(Q,b,c,s,x0,nMC,Ts);
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}
One could then plot the mean, standard deviation, a few Monte Carlo
runs, and an histogram at the final time using the following \matlab{}
commands
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5.25in}
\begin{lstlisting}
StMC=reshape(X(3,:,:),nMC,[])';   % extract only St values

% plot Mean, Std.Dev., and Monte Carlos runs
subplot(1,2,1)                    
h=plot(Ts,Xmean(3,:)','b-',...
       Ts,[Xmean(3,:)'-Xstd(3,:)',Xmean(3,:)'+Xstd(3,:)'],'b:',...
       Ts,StMC(:,1:10),'k-');
axis([0,Ts(end),0,1.1*x0(3)]);
legend(h([1,2,4]),'E[St]','E[St]\pm{}Std[St]','Monte Carlo runs')

% plot distribution
subplot(1,2,2)
[n,x]=hist(StMC(end,:),20);
barh(x,n)
axis([0,max(n),0,1.1*x0(3)]);
legend(sprintf('histogram at time %g',Ts(end)))
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}
These commands result in the plots in Figure~\ref{fig:MC}.
\begin{figure}
  \centering
  \includegraphics[width=.5\textwidth]{mc.eps}
  \caption{Results of the Monte Carlo simulations}
  \label{fig:MC}
\end{figure}

\medskip

The precise syntax of the \SDT{} commands above is explained in the
subsequent sections.

\cleardoublepage
\chapter{Moment Closure Methods}\label{se:closure}

Before describing the \SDT{} toolbox, we review moment closure and the
methods used in this software to achieve it.

\medskip

Consider a set of chemical species $X_1,X_2,\dots,X_n$ involved in a
set of chemical reactions and let us denote by $x\eqdef
(x_1,x_2,\dots,x_n)$ a vector containing their molecule counts. Given a
vector of integers $m\eqdef(m_1,m_2,\dots,m_n)$, we use the notation
$\mu^{(m)}$ to denote the following uncentered moment of $x$:
\begin{align*}
  \mu^{(m)}\eqdef\E[ x_1^{m_1} x_2^{m_2}\cdots x_n^{m_n}].
\end{align*}
Such moment if said to be of order $\sum_i m_i$. With $n$ species there
are exactly $n$ first order moments $\E[x_i]$, $\forall i
\in\{1,2,\dots,n\}$, which are just the means; $n(n-1)/2$ second order
moments $\E[x_i^2]$, $\forall i$ and $\E[x_i x_j]$, $\forall i\neq j$,
which can be used to compute variances and covariance; $n(n-1)(n-2)/6$
third order moments; and so on.

\medskip

It was show in \cite{SinghHespanhaDec06}, that if we construct a vector
$\mu$ containing all the uncentered moments of $x$ up to some order
$k$, the evolution of $\mu$ is determined by a differential equation of
the form
\begin{align}\label{eq:exact}
  \dot\mu &= A \mu+B \bar\mu, \quad \mu\in\R^K,\;\bar \mu\in\R^{\bar K}
\end{align}
where $A$ and $B$ are appropriately defined matrices and $\bar\mu$ is a
vector containing moments of order larger\footnote{When one does not
  include in $\mu$ the zero-order moment $\mu^{(0)}=1$, this term
  will appear in $\bar\mu$.} than $k$. The equation
\eqref{eq:exact} is exact and we call it the \emph{(exact) $k$-order
  moment dynamics} and the integer $k$ is called the \emph{order of
  truncation}.  Note that the dimension $K$ of \eqref{eq:exact} is
always larger than $k$ since there are many moments of each order. In
fact, in general $K$ is of order $n^k$.

\medskip

When all chemical reactions have only one reactant, the term $B\bar\mu$
does not appear in \eqref{eq:exact} and we say that the exact moment
dynamics are \emph{closed}. However, when at least one chemical
reaction has 2 or more reactants, then the term $B\bar\mu$ appears and
we say that the moment dynamics are \emph{open} since \eqref{eq:exact}
depends on the moments in $\bar\mu$, which are not part of the state
$\mu$.  When all chemical reactions are elementary (i.e., with at most
2 reactants), then all moments in $\bar\mu$ are exactly of order $k+1$.

\medskip

\emph{Moment closure} is a procedure by which one approximates the
exact (but open) moment dynamics \eqref{eq:exact} by an approximate
(but now closed) equation of the form
\begin{align}\label{eq:approx-1}
  \dot\nu = A \nu+B \varphi(\nu),\quad \nu\in\R^K
\end{align}
where $\varphi(\nu)$ is a column vector that approximates the moments in
$\bar\mu$. The function $\varphi(\nu)$ is called the moment closure function
and \eqref{eq:approx-1} is called the \emph{approximate $k$th order
  moment dynamics}.  The goal of any moment closure method is to
construct $\varphi(\nu)$ so that the solution $\nu$ to \eqref{eq:approx-1} is
close to the solution $\mu$ to \eqref{eq:exact}.

\medskip

Some moment closure methods approximate the exact moments dynamics
\eqref{eq:exact} by a closed equation of larger order, such as in
\begin{subequations}\label{eq:approx-2}
  \begin{align}\label{eq:approx-2-phi}
    \dot\phi &= \psi(\phi), \quad \phi\in\R^N, \\
    \label{eq:approx-2-nu}
    \dot\nu &= A \nu+B \varphi(\phi,\nu), \quad \nu\in\R^K,
  \end{align}
\end{subequations}
where one now approximates $\bar\mu$ by the function
$\varphi(\phi,\nu)$ that is allowed to also depend on the state $\phi$
of an additional dynamic system. Often $\varphi(\phi,\nu)$ can be made
linear in $\nu$. In this case, once $\phi$ reaches a steady state, the
$\nu$ dynamics became linear and time-invariant.

\medskip

There are three main approaches to construct the moment closure
function $\varphi(\cdot)$:
\begin{enumerate}
\item \emph{Matching-based methods} directly attempt to match the
  solutions to \eqref{eq:exact} and
  \eqref{eq:approx-1} [or \eqref{eq:approx-2}].

\item \emph{Distribution-based methods} construct $\varphi(\cdot)$ by
  making ``reasonable'' assumptions on the statistical distribution of
  the molecule counts vector $x$.

\item \emph{Large volume methods} construct $\varphi(\cdot)$ by assuming that
  reactions take place on a large volume.
\end{enumerate}
It is important to emphasize that this classification is about methods
to \emph{construct} moment closure. It turns out that sometimes
different methods lead to the same moment closure function $\varphi(\cdot)$.

\medskip

\begin{table*}[t!]
  \caption{Which moment closure to use?\label{ta:which}}
  \centering\codesize
  \begin{tabular}{||l||p{1.8in}|p{2.0in}|p{1.2in}||}\hline\hline
     & distributions have low variability (i.e., low standard deviations when
    compared to the mean) or are fairly symmetric & 
    distributions have large standard deviations when compared to the
    mean, but populations do not become zero with high probability &
    populations can become zero with high probability\\\hline\hline
    accuracy& zero cumulant closure &derivative matching
    closure& no good solution (yet) \\\hline
    simple dynamics& quasi-deterministic closure or Van Kampen's linear noise approximation & derivative
    matching closure (but will not be very simple)& no good solution (yet) \\\hline\hline
  \end{tabular}
\end{table*}
In the remainder of this section we discuss several methods to
construct the moment closure function. We shall see that the choice of which method to use
depends on the type of system (e.g., how population means compare with
standard deviations for the system considered) and also on the primary
goal in constructing the approximate moment dynamics (e.g., how
important is accuracy versus simplicity of the equations).
Table~\ref{ta:which} summarizes some rules of thumb on the choice of
which approximation to use.

\section{Derivative Matching}\label{se:derivative-matching}

\emph{Derivative matching} is a matching-based method for moment
closure described in \cite{SinghHespanhaDec06}. It uses moment closure
functions $\varphi(\cdot)$ in \eqref{eq:approx-1} whose entries are
\emph{separable}, i.e., of the form
\begin{align*}
  \nu_1^{\gamma_1} \nu_2^{\gamma_2} \cdots\nu_n^{\gamma_n}.
\end{align*}
The coefficients $\gamma_i\in\R$ are then computed to make the relative
error
\begin{align*}
  \frac{\Big\|\Deriv{^\ell\nu}{t^\ell}-\Deriv{^\ell\mu}{t^\ell}\Big\|}{\Big\|\Deriv{^\ell\mu}{t^\ell}\Big\|}
\end{align*}
as small as possible for molecule counts larger than one. Somewhat
surprisingly, this minimization leads to explicit formulas for the
moment closure functions $\varphi(\cdot)$ that do not depend on the
reaction parameters \cite{SinghHespanhaDec06}.

\section{Zero Cumulants}\label{se:zero-cumulants}

\emph{Zero cumulants} is a distribution-based method for moment
closure that finds the $k$th order moment closure function $\varphi(\cdot)$ in
\eqref{eq:approx-1} by assuming that all multi-variable cumulants of
the population $x$ with order larger than $k$ are negligible. This
makes the distribution of $x$ ``as close as possible'' to a
multi-variable Gaussian distribution, which has all cumulants or order
higher than two equal to zero.

\medskip

To construct zero cumulant closures, one uses the fact that the
cumulant $\kappa^{(m)}$ can be expressed as
\begin{align}\label{eq:kappa-mu}
  \kappa^{(m)}= \mu^{(m)}+\sum_{\sum \bar m_i<\sum m_i}  \alpha_{\bar m}\, \mu^{(\bar m)},
\end{align}
where the summation is over moments $\mu^{(\bar m)}$ of order strictly
smaller than $\sum m_i$ and the $\alpha_{\bar m}$ are appropriately selected
constants (c.f.~Appendix~\ref{se:cumulants}). This shows that the
cumulant $\kappa^{(m)}$ depends only on the moment $\mu^{(m)}$ and
lower-order moments $\mu^{(\bar m)}$, so by setting $\kappa^{(m)}=0$ one
obtains an expression for $\mu^{(m)}$ as a function of lower-order
moments.

\medskip

The procedure to compute the \emph{zero-cumulants} moment closure function
$\varphi(\cdot)$ consists of setting to zero all cumulants corresponding to
the moments that do not appear in $\mu$ and then solving the equations
\eqref{eq:kappa-mu} for the moments in $\bar\mu$.

\section{Low Dispersion}\label{se:low-dispersion}

\emph{Low dispersion} is a distribution-based method for moment
closure that finds the moment closure function $\varphi(\cdot)$ in
\eqref{eq:approx-1} by assuming that the distributions of the
populations are tightly clustered around their means, with standard
deviations much smaller than the means.  Specifically, for the $k$th
order moment closure one assumes that the normalized centered moments
of order larger than $k$ are much smaller than one.  We recall that
given a vector of integers $m\eqdef(m_1,m_2,\dots,m_n)$, the
corresponding \emph{normalized centered moment} is defined by
\begin{align*}
  \eta^{(m)}\eqdef\E\Big[\Big(\frac{x_1-\E[x_1]}{\E[x_1]}\Big)^{m_1}
  \Big(\frac{x_2-\E[x_2]}{\E[x_2]}\Big)^{m_2} \cdots
  \Big(\frac{x_n-\E[x_n]}{\E[x_n]}\Big)^{m_n}
  \Big].
\end{align*}
Such moment if said to be of order $\sum_i m_i$.
%
%\medskip
%
For fairly symmetric distributions the odd-order moments can be quite
small and therefore this technique is especially useful for even-order
moment closures for which the odd-order higher moments can be safely
neglected.

\medskip

To construct low dispersion closures, one uses the fact that an
uncentered moment $\mu^{(m)}$ can be expressed in terms of the
normalized centered moment as follows
\begin{align*}
  \mu^{(m)} &=\E[x_1]^{m_1}\E[x_2]^{m_2}\cdots\E[x_n]^{m_n}
 \Big(1+\eta^{(m)}+\sum_{2\leq \sum\bar m_i< \sum m_i}\beta_{\bar m}\,\eta^{(\bar m)}\Big),
\end{align*}
where the summation is over moments $\eta^{(\bar m)}$ of order two or
larger and strictly smaller than $\sum m_i$, and the $\beta_{\bar m}$
are appropriately selected nonnegative constants
(cf.~Appendix~\ref{se:centered}). When a particular normalized
centered moment $\eta^{(m)}$ is much smaller than one, we have that
\begin{align}\label{eq:mu-approx-low}
  \mu^{(m)} \approx \E[x_1]^{m_1}\E[x_2]^{m_2}\cdots\E[x_n]^{m_n}
 \Big(1+\sum_{2\leq \sum\bar m_i< \sum m_i}\beta_{\bar m}\,\eta^{(\bar m)}\Big),
\end{align}
which allows one to express the uncentered moment $\mu^{(m)}$ solely in
terms of normalized centered moments $\eta^{(\bar m)}$ of order strictly
smaller than $\sum m_i$.  On the other hand, we can express all these
normalized centered moments as linear combinations of the uncentered
moments of order strictly smaller than $\sum m_i$ as follows
\begin{align}\label{eq:eta-mu-low}
  \eta^{(\bar m)}
  &=\frac{\mu^{(\bar m)}}{\E[x_1]^{\bar m_1}\E[x_2]^{\bar m_2}\cdots\E[x_n]^{\bar m_n}}
  + \sum_{\sum\tilde m_i< \sum m_i}\gamma_{\tilde m}\,\frac{\mu^{(\tilde m)}}{\E[x_1]^{\tilde m_1}\E[x_2]^{\tilde m_2}\cdots\E[x_n]^{\tilde m_n}},
\end{align}
where the summation is over uncentered moments $\mu^{(\tilde m)}$ of
order strictly smaller than $\sum \bar m_i$ and the $\gamma_{\tilde m}$
are appropriately selected constants (cf.~Appendix~\ref{se:centered}).

\medskip

The procedure to compute the \emph{low dispersions} moment closure
function $\varphi(\cdot)$ in \eqref{eq:approx-1} thus consists of using
\eqref{eq:mu-approx-low} and \eqref{eq:eta-mu-low} to approximate any
moment that does not appear in $\mu$ as a linear combination of the
moments in $\mu$. Note however that the coefficients of these linear
combinations will depend on monomials of the form $\E[x_1]^{\hat
  m_1}\E[x_2]^{\hat m_2}\cdots\E[x_n]^{\hat m_n}$, with all the $\hat
m_i\geq 0$ and therefore the moment closure function $\phi(\cdot)$ will be
polynomial but nonlinear on $\mu$.

\paragraph{Relationship with zero-cumulants closure} For second order
moment closure ($k=2$) one sets to zero 3th-order normalized centered
moments, which is equivalent to setting to zero the 3th-order
cumulants. Therefore for 2nd-order closures, zero cumulant and low
dispersion coincide.

\section{Quasi-Deterministic}\label{se:quasi-deterministic}

\emph{Quasi-deterministic} is a distribution-based method for moment
closure that finds the moment closure function $\varphi(\cdot)$ in
\eqref{eq:approx-2} by assuming that the distributions of the
populations are tightly clustered around the solution $\phi$ to the
deterministic dynamics
\begin{align}\label{eq:qd-phi}
  \dot \phi&=A_\mrm{det} \phi+B_\mrm{det} \psi(\phi), & \phi\eqdef(\phi_1,\phi_2,\dots,\phi_n)\in\R^n,
\end{align}
which are obtained by assuming that each $\phi_i\eqdef x_i$ is
deterministic and therefore
\begin{align*}
  \E[\phi_i \phi_j] = \E[\phi_i] \E[\phi_j] = \phi_i \phi_j.
\end{align*}
Specifically, for the $k$th order moment closure one assumes that the
quasi-deterministic normalized centered moments of order larger than
$k$ are much smaller than one.  Given a vector of integers
$m\eqdef(m_1,m_2,\dots,m_n)$, the corresponding
\emph{quasi-deterministic normalized centered moment} is defined by
\begin{align*}
  \hat \eta^{(m)}\eqdef\E\Big[\Big(\frac{x_1-\phi_1}{\phi_1}\Big)^{m_1}
  \Big(\frac{x_2-\phi_2}{\phi_2}\Big)^{m_2} \cdots
  \Big(\frac{x_n-\phi_n}{\phi_n}\Big)^{m_n}
  \Big].
\end{align*}
Such moment if said to be of order $\sum_i m_i$.

\medskip

To construct quasi-deterministic closures, one uses the fact that an
uncentered moment $\mu^{(m)}$ can be expressed in terms of the
quasi-deterministic normalized centered moment as follows
\begin{align*}
  \mu^{(m)} &=\phi_1^{m_1}\phi_2^{m_2}\cdots\phi_n^{m_n}
 \Big(1+\hat\eta^{(m)}+\sum_{1\leq \sum\bar m_i< \sum m_i}\beta_{\bar m}\,\hat\eta^{(\bar m)}\Big),
\end{align*}
where the summation is over moments $\eta^{(\bar m)}$ of order one or
larger and strictly smaller than $\sum m_i$, and the $\beta_{\bar m}$
are appropriately selected nonnegative constants
(cf.~Appendix~\ref{se:ap-quasi-deterministic}). When a particular
quasi-deterministic normalized centered moment $\hat\eta^{(m)}$ is
much smaller than one, we have that
\begin{align}\label{eq:mu-approx-qd}
  \mu^{(m)} \approx \phi_1^{m_1}\phi_2^{m_2}\cdots\phi_n^{m_n}
 \Big(1+\sum_{2\leq \sum\bar m_i< \sum m_i}\beta_{\bar m}\,\hat\eta^{(\bar m)}\Big),
\end{align}
which allows one to express the uncentered moment $\mu^{(m)}$ solely in
terms of quasi-deterministic normalized centered moments $\hat
\eta^{(\bar m)}$ of order strictly smaller than $\sum m_i$.  On the other
hand, we can express all these quasi-deterministic normalized centered
moments as linear combinations of the uncentered moments of order
strictly smaller than $\sum m_i$ as follows
\begin{align}\label{eq:eta-mu-qd}
  \hat \eta^{(\bar m)}
  &=\frac{\mu^{(\bar m)}}{\phi_1^{\bar m_1}\phi_2^{\bar m_2}\cdots\phi_n^{\bar m_n}}
  + \sum_{\sum\tilde m_i< \sum m_i}\gamma_{\tilde m}\,\frac{\mu^{(\tilde m)}}{\phi_1^{\tilde m_1}\phi_2^{\tilde m_2}\cdots\phi_n^{\tilde m_n}},
\end{align}
where the summation is over uncentered moments $\mu^{(\tilde m)}$ of
order strictly smaller than $\sum \bar m_i$ and the $\gamma_{\tilde m}$
are appropriately selected constants
(cf.~Appendix~\ref{se:ap-quasi-deterministic}).

\medskip

The procedure to compute the \emph{quasi-deterministic} moment
closure function $\varphi(\cdot)$ in \eqref{eq:approx-2} thus consists of using
\eqref{eq:mu-approx-low} and \eqref{eq:eta-mu-low} to approximate any
moment that does not appear in $\mu$ as a linear combination of the
moments in $\mu$. The coefficients of these linear combinations will
depend on monomials of the form $\phi_1^{\hat m_1}\phi_2^{\hat
  m_2}\cdots\phi_n^{\hat m_n}$, with all the $\hat m_i\geq 0$ and
therefore the moment closure function will be linear on $\mu$ for a
fixed $\phi$. This means that the approximate
dynamics in~\eqref{eq:approx-2} are of the form
\begin{align}\label{eq:approx-3}
    \dot \phi&=A_\mrm{det} \phi+B_\mrm{det} \psi(\phi), \quad \phi\in\R^n, &
    \dot\nu &= A_\nu(\phi)\nu+c_\nu(\phi), \quad \nu\in\R^K,
\end{align}
and, when $\phi$ reaches a steady state value, the $\nu$ dynamics become
linear.

\paragraph{Relationship with low dispersion closure} In general the
normalized centered moment are smaller than their quasi-deterministic
version and therefore whenever quasi-deterministic moment closure
provides a good approximation, one should expect low-dispersion moment
closure to do at least as well. However, quasi-deterministic moment
closure has the advantage that it results in moment dynamics that are
``almost'' linear and therefore generally easier to analyze.

\section{Van Kampen's Linear Noise Approximation}\label{se:lna}

Van Kampen's \emph{Linear Noise Approximation} is developed in
\cite[Chapter X]{Kampen01} and can be applied when the matrices
$A,B$ in \eqref{eq:exact} depend on some parameter $V$ that can be
assumed large, i.e., when we have
\begin{align*}
  \dot\mu &= A(V) \mu+B(V) \bar\mu, \quad \mu\in\R^K,
\end{align*}
with $V$ large.  This form of moment closure results in a system of
the form \eqref{eq:approx-2} and is exact in the limit as $V\to\infty$.
Typically, $V$ is the volume on which the chemical reactions take
place.

\medskip

To construct \eqref{eq:approx-2}, one starts by choosing $\bar\phi$ to
satisfy the \emph{deterministic large-volume dynamics}
\begin{align}\label{eq:lna-phi}
  \dot{\bar\phi}&=A_\mrm{det} \bar\phi+B_\mrm{det} \psi(\bar\phi), & \bar\phi\eqdef(\bar\phi_1,\bar\phi_2,\dots,\bar\phi_n)\in\R^n
\end{align}
which are obtained by assuming that each $\bar\phi_i\eqdef x_i/V$ is deterministic
and therefore
\begin{align*}
  \E[\bar\phi_i \bar\phi_j] = \E[\bar\phi_i] \E[\bar\phi_j] = \bar\phi_i \bar\phi_j
\end{align*}
and also by making $V\to\infty$.

\medskip

Regarding the vector $\bar\phi$ in \eqref{eq:lna-phi} as a deterministic
approximation to the stochastic vector $x/V$, motivates defining the
following stochastic perturbation vector
$\chi\eqdef(\chi_1,\chi_2,\dots,\chi_n)$, with
\begin{align}\label{eq:chi-x}
  \chi_i &\eqdef \frac{x_i-V \bar\phi_i}{V^\frac{1}{2}}
  \eqv x_i = V \bar\phi_i + V^\frac{1}{2} \chi_i,
\end{align}
where the normalization by $V^\frac{1}{2}$ will be needed to keep the
moments of $\chi$ bounded as $V\to\infty$. Given a vector of integers
$m\eqdef(m_1,m_2,\dots,m_n)$, we use the notation $\xi^{(m)}$ to
denote the following uncentered moment of $\chi$:
\begin{align*}
  \xi^{(m)}\eqdef\E[ \chi_1^{m_1} \chi_2^{m_2}\cdots \chi_n^{m_n}].
\end{align*}
The moments of $x$ and $\chi$ are related by
\begin{align}\label{eq:mu-vs-xi}
  \mu^{(m)}%=\E[x_1^{m_1} x_2^{m_2}\cdots x_n^{m_n}]
  = \E[(V \bar\phi_1 + V^\frac{1}{2} \chi_1)^{m_1} \cdots (V \bar\phi_n + V^\frac{1}{2} \chi_n)^{m_n}]
  = V^{\sum_i m_i}\sum_{\sum\bar m_i\leq \sum m_i} \frac{\alpha_{\bar m}}{V^{\frac{\bar m}{2}}} \xi^{(\bar m)}
  = V^{\sum_i m_i}\Big(\bar\phi^{(m)} + \cdots +\frac{1}{V^{\frac{\sum_i m_i}{2}}}\xi^{(m)}\Big),
\end{align}
where the summation is over moments $\xi^{(\bar m)}$ of order up to
$\sum m_i$ and the $\alpha_{\bar m}$ are appropriately selected
constants. 

\medskip

Computing the (exact) moment dynamics for $\xi$, one
obtains\footnote{Assuming that $\xi$ remain bounded as $V\to\infty$, for any
  moment $\mu^{(m)}$ in $\bar\mu$, we have that for sufficiently large
  $V$,
  \begin{align*}
    \mu^{(m)} \approx V^{\sum_i m_i}\sum_{\bar m\leq k} \frac{\alpha_{\bar m}}{V^{\frac{\bar m}{2}}} \xi^{(\bar m)}
  \end{align*}
  where we only kept in the summation the terms order up to $k$ that
  are in $\mu$. Since all the $\xi^{(\bar m)}$, $\bar m\leq k$ can be
  expressed as linear combinations of the moments in $\mu$, we have
  obtained an expression for the moments in $\bar\mu $ as a linear
  combination of the moments in $\mu$ (generally with time varying
  coefficient that depend on $\bar\phi$). This approach only make
  $V\to\infty$ where this is absolutely needed and therefore could lead to
  better results than taking the limit in \eqref{eq:xi-2}. However, in
  the current \SDT{} implementation we followed closely \cite[Chapter
  X]{Kampen01} and took the limit in \eqref{eq:xi-2}.}
\begin{align}\label{eq:xi-1}
  \dot\xi &= A_\xi(V,\bar\phi) \xi+B_\xi(V,\bar\phi) \bar\xi, \quad \xi\in\R^K,
\end{align}
where $\xi$ and $\bar\xi$ contain the moments of $\chi$ corresponding to
the moments of $x$ in $\mu$ and $\bar\mu$, respectively.  For
elementary reactions with reaction rates that depend on the volume as
follows:
\begin{align*}
  \emptyset  & \xrightarrow{\text{rate } = c V} * &
  X  & \xrightarrow{\text{rate } = c X} *  & 
  2 X  & \xrightarrow{\text{rate } = c\frac{X(X-1)}{V}} *  &
  X+Y  & \xrightarrow{\text{rate } = c\frac{X Y}{V}} * 
\end{align*}
the open system \eqref{eq:xi-1} converges as $V\to\infty$ to a
closed system of the form
\begin{align}\label{eq:xi-2}
  \dot\xi &= A_\xi(V,\bar\phi) \xi+B_\xi(V,\bar\phi) \bar\xi 
  \xrightarrow{V\to\infty} A_\xi(\infty,\bar\phi) \xi, \quad \xi\in\R^K.
\end{align}
Since the moments $\mu$ and $\xi$ are related through
\eqref{eq:mu-vs-xi}, one can use \eqref{eq:xi-2} to obtain a closed
equation for $\mu$ as in \eqref{eq:approx-2}. Moreover, this
equation will be linear in $\mu$, leading to approximate dynamics
similar to \eqref{eq:approx-3}.

\medskip

One could have done the derivation above in a more constructive way
without pre-specifying the dynamics for $\bar\phi$ by
\eqref{eq:lna-phi}. In this case, since we have no expression to
replace for $\dot{\bar\phi}$, we would have obtained for the exact
dynamics for $\xi$ an expression of the form
\begin{align*}
  \dot\xi &= A_\xi(V,\bar\phi) \xi+B_\xi(V,\bar\phi)
  \bar\xi+C_\xi(V,\bar\phi)\dot{\bar\phi}, \quad \xi\in\R^K.
\end{align*}
If then tried to make $V\to\infty$, we would observe that to obtain a
finite right-hand side we would need $\dot{\bar\phi}$ to satisfy
precisely \eqref{eq:lna-phi}.

\paragraph{Relationship with quasi-deterministic closure} The
deterministic equations \eqref{eq:qd-phi} and \eqref{eq:lna-phi}
differ by two facts:
\begin{enumerate}
\item the state in \eqref{eq:lna-phi} was normalized through a division by
  the volume,
\item in \eqref{eq:lna-phi} we took the limit as $V\to\infty$.
\end{enumerate}
For elementary reactions with molecule counts much larger than one,
taking the limit as $V\to\infty$ has almost no effect and we essentially
have $\phi=V\bar \phi$. In this case, 
\begin{align*}
  \xi^{(m)} &= \E\Big[\Big(\frac{x_1-\phi_1}{V^\frac{1}{2}}\Big)^{m_1}
  \Big(\frac{x_2-\phi_2}{V^\frac{1}{2}}\Big)^{m_2} \cdots
  \Big(\frac{x_n-\phi_n}{V^\frac{1}{2}}\Big)^{m_n}
  \Big]
  =\frac{\phi_1^{m_1}\phi_2^{m_2}\cdots\phi_n^{m_n}}{V^{\frac{\sum_i m_i}{2}}}\hat \eta^{(m)}.
\end{align*}
So setting a quasi-deterministic moments to zero $\hat \eta^{(m)}$ is
equivalent to setting to zero the corresponding uncentered moment
$\xi^{(m)}$ of $\chi$. This means that we can view the
quasi-deterministic closure as taking the Van Kampen equations
\eqref{eq:lna-phi} and \eqref{eq:xi-2} and simply setting $\bar\xi$ in
\eqref{eq:xi-2} to zero, without ignoring other terms that would also
disappear as $V\to\infty$. Since we are keeping more terms of the exact
equations, with quasi-deterministic closure one often obtains
more accurate results then with Van Kampen's linear noise
approximation.

\section{Example}\label{se:example}

In this section we present the different moment closure dynamics for the
network of chemical reactions considered in \cite[p.~263]{Kampen01}.
This network of three elementary chemical reactions is described by
the following \texttt{.net} file:
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5.25in}
\begin{lstlisting}[language=net,columns={[l]flexible},basewidth={.6em,.45em}]
species:
        X stochastic;     % number of X molecules
        Y stochastic;     % number of Y molecules
parameters:
        V          = 20;  % volume
        phiA "\phi_A" = 5;   % concentration of A (fixed)
        al "\alpha"     = 10;
        be "\beta"     = 20;
        ga "\gamma"     = 30;
reactions:
        rate = al*phiA*V;      {X} > {X+1};          %   A -> X
        rate = ga*X*(X-1)/V;   {X,Y} > {X-2,Y+1};    % 2 X -> Y
        rate = be*Y;           {Y} > {Y-1};          %   Y -> B
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}  
The exact 2nd-order moment dynamics for this system are given by
\begin{align*}
  \Deriv{}{t}
  \smatt{\E[X]\\
    \E[Y]\\
    \E[X^2]\\
    \E[X Y]\\
    \E[Y^2]}=
  \smatt{2{\frac{\gamma}{V}}&0&-2{\frac{\gamma}{V}}&0&0\\
    -{\frac{\gamma}{V}}&-\beta&{\frac{\gamma}{V}}&0&0\\
    2 \alpha \phi_A V-4{\frac{\gamma}{V}}&0&8{\frac{\gamma}{V}}&0&0\\
    2{\frac{\gamma}{V}}&\alpha \phi_A V&-3{\frac{\gamma}{V}}&-\beta+2{\frac{\gamma}{V}}&0\\
    -{\frac{\gamma}{V}}&\beta&{\frac{\gamma}{V}}&-2{\frac{\gamma}{V}}&-2 \beta}
  \smatt{\E[X]\\
    \E[Y]\\
    \E[X^2]\\
    \E[X Y]\\
    \E[Y^2]}+
  \smatt{\alpha \phi_A V&0&0\\
    0&0&0\\
    \alpha \phi_A V&-4{\frac{\gamma}{V}}&0\\
    0&{\frac{\gamma}{V}}&-2{\frac{\gamma}{V}}\\
    0&0&2{\frac{\gamma}{V}}}
  \smatt{1\\
    \E[X^3]\\
    \E[X^2 Y]}.
\end{align*}
We now list the approximate 2nd-order moment dynamics for this system
obtained using the different methods:
\begin{enumerate}
\item Derivative matching:
  \begin{align*}
    \smatt{\E[X^3]\\
      \E[X^2 Y]}\approx
    \smatt{{\frac{{\E[X^2]}^{3}}{{\E[X]}^{3}}}\\
      {\frac{\E[X^2]{\E[X Y]}^{2}}{{\E[X]}^{2}\E[Y]}}}
  \end{align*}

\item Zero cumulants and low dispersion:
  \begin{align*}
    \smatt{\E[X^3]\\
      \E[X^2 Y]}\approx
    \smatt{3 \E[X^2]\E[X]-2{\E[X]}^{3}\\
      \E[X^2]\E[Y]+2 \E[X]\E[X Y]-2{\E[X]}^{2}\E[Y]}
  \end{align*}
  
\item  Quasi deterministic:
  \begin{align*}
    \Deriv{}{t}\smatt{\phi_{X}\\\phi_{Y}}
    &=\smatt{2{\frac{\gamma \phi_{X}}{V}}-2{\frac{\gamma{\phi_{X}}^{2}}{V}}+\alpha \phi_A V\\
      -{\frac{\gamma \phi_{X}}{V}}-\beta \phi_{Y}+{\frac{\gamma{\phi_{X}}^{2}}{V}}}\\
    \smatt{\E[X^3]\\
      \E[X^2 Y]}&\approx
    \smatt{{\phi_{X}}^{3}-3{\phi_{X}}^{2}\E[X]+3 \phi_{X}\E[X^2]\\
      {\phi_{X}}^{2}\phi_{Y}+\phi_{Y}\E[X^2]-2 \phi_{X}\phi_{Y}\E[X]-{\phi_{X}}^{2}\E[Y]+2 \phi_{X}\E[X Y]}
  \end{align*}
  leading to
  \begin{align*}
    \Deriv{}{t}\smatt{\E[X]\\
      \E[Y]\\
      \E[X^2]\\
      \E[X Y]\\
      \E[Y^2]}\approx
    \smatt{2{\frac{\gamma}{V}}&0&-2{\frac{\gamma}{V}}&0&0\\
      -{\frac{\gamma}{V}}&-\beta&{\frac{\gamma}{V}}&0&0\\
      2 \alpha \phi_A V-4{\frac{\gamma}{V}}+12{\frac{\gamma{\phi_{X}}^{2}}{V}}&0&8{\frac{\gamma}{V}}-12{\frac{\gamma \phi_{X}}{V}}&0&0\\
      2{\frac{\gamma}{V}}-3{\frac{\gamma{\phi_{X}}^{2}}{V}}+4{\frac{\gamma \phi_{X}\phi_{Y}}{V}}&\alpha \phi_A V+2{\frac{\gamma{\phi_{X}}^{2}}{V}}&-3{\frac{\gamma}{V}}+3{\frac{\gamma \phi_{X}}{V}}-2{\frac{\gamma \phi_{Y}}{V}}&-\beta+2{\frac{\gamma}{V}}-4{\frac{\gamma \phi_{X}}{V}}&0\\
      -{\frac{\gamma}{V}}-4{\frac{\gamma \phi_{X}\phi_{Y}}{V}}&\beta-2{\frac{\gamma{\phi_{X}}^{2}}{V}}&{\frac{\gamma}{V}}+2{\frac{\gamma \phi_{Y}}{V}}&-2{\frac{\gamma}{V}}+4{\frac{\gamma \phi_{X}}{V}}&-2 \beta}
    \smatt{\E[X]\\
      \E[Y]\\
      \E[X^2]\\
      \E[X Y]\\
      \E[Y^2]}
    +\smatt{\alpha \phi_A V\\
      0\\
      \alpha \phi_A V-4{\frac{\gamma{\phi_{X}}^{3}}{V}}\\
      {\frac{\gamma{\phi_{X}}^{3}}{V}}-2{\frac{\gamma{\phi_{X}}^{2}\phi_{Y}}{V}}\\
      2{\frac{\gamma{\phi_{X}}^{2}\phi_{Y}}{V}}}
  \end{align*}

\item  Van Kampen's linear noise approximation:
  \begin{align*}
    \Deriv{}{t}\smatt{\phi_{X}\\\phi_{Y}}
    &=\smatt{-2 \gamma{\phi_{X}}^{2}+\alpha \phi_A\\
-\beta \phi_{Y}+\gamma{\phi_{X}}^{2}}
    \\
    \Deriv{}{t}\smatt{\xi_{X}\\
      \xi_{Y}\\
      \xi_{X^2}\\
      \xi_{X Y}\\
      \xi_{Y^2}}
    &\approx
    \smatt{-4 \gamma \phi_{X}&0&0&0&0\\
      2 \gamma \phi_{X}&-\beta&0&0&0\\
      0&0&-8 \gamma \phi_{X}&0&0\\
      0&0&2 \gamma \phi_{X}&-\beta-4 \gamma \phi_{X}&0\\
      0&0&0&4 \gamma \phi_{X}&-2 \beta}
    \smatt{\xi_{X}\\
      \xi_{Y}\\
      \xi_{X^2}\\
      \xi_{X Y}\\
      \xi_{Y^2}}
    +\smatt{0\\
      0\\
      4 \gamma{\phi_{X}}^{2}+\alpha \phi_A\\
      -2 \gamma{\phi_{X}}^{2}\\
      \gamma{\phi_{X}}^{2}+\beta \phi_{Y}}
  \end{align*}
  leading to
  \begin{align*}
    \Deriv{}{t}\smatt{\E[X]\\
      \E[Y]\\
      \E[X^2]\\
      \E[X Y]\\
      \E[Y^2]}\approx
    \smatt{-4 \gamma \phi_{X}&0&0&0&0\\
      2 \gamma \phi_{X}&-\beta&0&0&0\\
      4 \gamma V{\phi_{X}}^{2}+2 \alpha \phi_A V&0&-8 \gamma \phi_{X}&0&0\\
      -\gamma V{\phi_{X}}^{2}&2 \gamma V{\phi_{X}}^{2}+\alpha \phi_A V&2 \gamma \phi_{X}&-\beta-4 \gamma \phi_{X}&0\\
      0&-2 \gamma V{\phi_{X}}^{2}&0&4 \gamma \phi_{X}&-2 \beta}
    \smatt{\E[X]\\
      \E[Y]\\
      \E[X^2]\\
      \E[X Y]\\
      \E[Y^2]}
    +\smatt{2 \gamma V{\phi_{X}}^{2}+\alpha \phi_A V\\
      -\gamma V{\phi_{X}}^{2}\\
      4 \gamma V{\phi_{X}}^{2}+\alpha \phi_A V\\
      -2 \gamma V{\phi_{X}}^{2}\\
      \gamma V{\phi_{X}}^{2}+\beta V \phi_{Y}}
  \end{align*}
\end{enumerate}
Figure~\ref{fig:VanKampen} compares the accuracy of the different
moment closure methods for a low volume ($V=2$) and a high volume
($V=20$).  For the larger volume all moment closure techniques provide
a very good match with Monte Carlo results, but for the smaller volume
derivative matching produces the most accurate results even with only
a second order truncation. These results are fairly typical.

\begin{figure}
  \centering
  \subfigure[volume $V=2$]{\includegraphics[height=.4\textheight]{VanKampen_V_2}}\\
  \subfigure[volume $V=20$]{\includegraphics[height=.4\textheight]{VanKampen_V_20}}\\
  
  \caption{Comparison of accuracy between different moment closure
    methods for the example described in Section~\ref{se:example} with
    two different volumes.  The legends show (i) the values of the
    mean $\pm$ one standard deviation at the final time (ii) a
    two-character string indicating the moment closure method, and
    (iii) an integer indicating the order of truncation.  The
    distributions, means and standard deviations in the right-most
    plots were obtained using 20,000 Monte Carlo simulations produced
    by \cite{StochKit}. The left-most plots include a typical Monte
    Carlo run.}
  \label{fig:VanKampen}
\end{figure}


\chapter{Specifying networks of chemical reactions}

In \SDT{}, networks of chemical reactions are specified by
``\texttt{.net} files'' whose syntax is described in this section.
These files are read using the \SDT{} function \lstinline{readNet()},
which is also described here.

{\lstset{language=net}

\section{The \texttt{.net} file}


A \texttt{.net} file can contain up to five sections:
\lstinline{species}, \lstinline{parameters},
\lstinline{substitutions}, \lstinline{reactions}, and
\lstinline{derivatives}. The \lstinline{substitutions} and
\lstinline{derivatives} sections are optional. Examples of
\texttt{.net} files were provided in Section~\ref{se:quick}. 


\paragraph{Species section}

The \lstinline{species} section describes the chemical species involved
in the network. It also specifies some assumptions that can be made in
analyzing the system. Following a \lstinline{species:} header, each
line is of the form:
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5in}
\begin{lstlisting}
species_name "species_latex_name" species_type initial_value;
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}
where
\begin{itemize}
\item \lstinline{species_name} stands for the symbol that identifies
  a chemical species. 
  
  The symbol should contain no spaces and be a valid name for a
  \matlab{} variable\footnote{\label{fo:MuPAD} The following names
    should be avoided as they are predefined MuPAD symbols:
    \lstinline{beta}, \lstinline{gamma}, \lstinline{psi},
    \lstinline{theta}, \lstinline{zeta}, \lstinline{D},
    \lstinline{E}, \lstinline{I}, \lstinline{O},
    \lstinline{Ei}, \lstinline{Ci}, \lstinline{Si}.}

\item \lstinline{species_latex_name} stands for an optional \LaTeX{}
  string that represents the species population to be used when
  expressions are converted to a \LaTeX{} form.

  This string appears quoted by " and should not include any other quotes.
  
\item \lstinline{species_type} stands for the assumed type for the
  species.  
  
  The species type specifies assumptions that can be made in analyzing
  the system. It can be one of the following keywords:
  \begin{itemize}
  \item \lstinline{constant} should be used when the population of
    the species remains constant;

  \item \lstinline{boolean} should be used when only 0 or 1 molecules
    can be present (typically genes that can either be active or
    inactive);

  \item \lstinline{deterministic} should be used when the expected number
    of molecules is much larger than its standard deviation and
    therefore stochastic effects can be neglected for this species;
    
  \item \lstinline{stochastic} should be used when no assumptions can
    be made about this species.
  \end{itemize}
  
  \textbf{Attention!} In principle one could declare all populations
  as \lstinline{stochastic}, but there are several advantages to
  providing more information about the populations:
  \begin{enumerate}
  \item Assuming that the assumptions are indeed true, the
    moment dynamics obtained will be more accurate.

  \item The order of the set of ODEs will generally be lower.
    
  \item For populations with low stochasticity (i.e., for which
    standard deviations are much smaller than the means), if one does
    not declare the population as \lstinline{deterministic} one runs the
    risk of obtaining slightly negative variances due to the moment
    closure approximation errors.
  \end{enumerate}
  
\item \lstinline{initial_value} is an optional constant that specifies the
  initial number of molecules for numerical simulations.
  
  \textbf{Attention!} Starting with initial conditions for which some
  stochastic species have exactly zero molecules can be problematic
  when using moment closure based on derivative matching (leading to
  divisions by zero or errors in the ODE solver). To avoid this, one
  can initialize the system with a small but positive population (say
  .01 molecules).
\end{itemize}

\paragraph{Parameters section}

The \lstinline{parameters} section declares all parameters that appear in
the \lstinline{reactions} and \lstinline{derivatives} sections and provides
default numerical values for the parameters .  Following a
\lstinline{parameters:} header, each line is of the form:
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5in}
\begin{lstlisting}
parameter_name "parameter_latex_name" = default_value;
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}
where
\begin{itemize}
\item \lstinline{parameter_name} stands for the symbol that identifies
  the parameter.
  
  The symbol should contain no spaces and be a valid name for a
  \matlab{} variable\footnote{See footnote \ref{fo:MuPAD}}.
  
\item \lstinline{parameter_latex_name} stands for an optional \LaTeX{}
  string that represents the parameter to be used when expressions are
  converted to a \LaTeX{} form.

  This string appears quoted by " and should not include any other quotes.
  
\item \lstinline{default_value} is a symbolic expression that stands for the
  default value of the parameter.
 
  All symbolic computations ignore this value and treat parameters as
  symbolic variables. Default values are only used when numerical
  values are needed.
\end{itemize}

\paragraph{Substitutions section}

The \lstinline{substitutions} section provides replacement rules that
should be used to simplify computations that involve the populations
of the different species.  Following a \lstinline{rules:}
header, each line is of the form:
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5in}
\begin{lstlisting}
{old_expression} > {new_expression};
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}
where
\begin{itemize}
\item \lstinline{old_expression} stands for a valid \matlab{} expression that
  should be replaced by \lstinline{new_expression}. The rule is
  applied multiple times until a ``fixed point'' is achieved.

  \medskip
  
  For example, if a species \lstinline{X} is known to have either 1 or 2
  individuals, then one should use the rule
  \begin{lstlisting}
    {X^2} > {3*X-2}
  \end{lstlisting}
  On the other hand, if the population of \lstinline{X} is known to be
  in the set $\{0,1,2\}$, then one should use
  \begin{lstlisting}
    {X^3} > {3*X^2-2*X}
  \end{lstlisting}
\end{itemize}

Currently, these rules are used in the following situations
\begin{itemize}
\item To simplify the computation of moment dynamics (prior to any
  moment closure) in the functions \lstinline{closureDynamics()} and
  \lstinline{momentDynamics()}.
  
\item To simplify the computation of expected values, standard
  deviations, and distributions in the functions
  \lstinline{getCMoments()}, \lstinline{plotCMoments}, and
  \lstinline{getDistribution()}.
\end{itemize}

These rules are applied \emph{before} processing any simplification
rules that arise from species being declared as \lstinline{boolean} or
\lstinline{deterministic}.

\paragraph{Reactions section}


The \lstinline{reactions} section describes the chemical reactions
involved in the network, including their stoichiometry and rates.
Following a \lstinline{reactions:} header, each line is of the form:
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5in}
\begin{lstlisting}
rate = rate_expr;  {list_species} > {post_reaction_counts};
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}
where
\begin{itemize}
\item \lstinline{rate_expr} stands for an expression describing the
  rate at which the reaction occurs. 
  
  The expression should be a \emph{polynomial}, possibly dependent on
  symbolic variables for which numerical values will be provided
  later. It should be a valid \matlab{} expression.
  
\item \lstinline{list_species} is a comma-separated list with the
  symbols of the chemical species whose stoichiometry changes in the
  reaction.  
  
\item \lstinline{post_reaction_counts} is a comma-separated list of
  expressions that specifies how the molecule counts for each symbol
  changes when the reaction takes place.  
\end{itemize}

\paragraph{Derivatives section}


The \lstinline{derivatives} section describes equations for possible
continuous evolutions for some/all of the species populations.
Following a \lstinline{derivatives:} header, each line is of the
form:
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5in}
\begin{lstlisting}
d/dt species_name += derivative_expr;
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}
where
\begin{itemize}
\item \lstinline{species_name} stands for the symbol of the chemical species
  whose derivative is provided.

\item \lstinline{derivative_expr} stands for an expression describing
  the rate of change of the population of \lstinline{species}. This
  change in the population is to be \emph{added} to the discrete rules
  specified in the \lstinline{reactions} section and to other rules for
  the same species that may be specified in the \lstinline{derivatives}
  section.

  The expression should be a \emph{polynomial}, possibly dependent on
  symbolic variables for which numerical values will be provided
  later. It should be a valid \matlab{} expression.
\end{itemize}

\paragraph{Comments}
Comments can be inserted anywhere in a \texttt{.net} file with the
prefix \lstinline{%}.
  
}

\section{Reading a \texttt{.net} file}

A \texttt{.net} file can be read using the following \SDT{} function:

\begin{command}{readNet()}
\begin{lstlisting}
net=readNet(filename)
\end{lstlisting}
\end{command}

\noindent
This \SDT{} function reads a \texttt{.net} file and stores its contents in
the \lstinline{net} structure.

\paragraph{Inputs and outputs:} 
\begin{itemize}
\item \lstinline{filename} specifies the name of the \texttt{.net} file
  to read
  
\item \lstinline{net} is a structure that describes the network in a
  format that is recognized by the \SDT{} toolbox. This structure has
  the following main fields:
  \begin{itemize}
  \item \lstinline{net.species} describes the species involved (from the
    \lstinline{species} section of the \texttt{.net} file);
  \item \lstinline{net.parameter} contains default numerical values for
    parameters that appear in the symbolic expressions in
    \lstinline{net.reactions} and \lstinline{net.raterule} (from the
    \lstinline{parameters} section of the \texttt{.net} file);
  \item \lstinline{net.substitutionrule} describes the replacement rules
    that should be used to simplify computations that involve
    populations (from the \lstinline{rules} section of the
    \texttt{.net} file).
  \item \lstinline{net.reaction} describes the chemical reactions (from
    the \lstinline{reactions} section of the \texttt{.net} file).
  \item \lstinline{net.raterule} describes the continuous rates of
    change for the populations (from the \lstinline{derivatives} section
    of the \texttt{.net} file).
  \end{itemize}
\end{itemize}

\chapter{Computing moment dynamics}

In \SDT{}, the exact and approximate moment dynamics for a network of
chemical reactions are characterized by the \lstinline{mdyn} structure
described in this section.
 
\section{Creating a \lstinline{mdyn} structure}

The approximate moment dynamics described in Section~\ref{se:closure}
can be computed using the \SDT{} function
\lstinline{closureDynamics()}.

\begin{command}{closureDynamics()}
\begin{lstlisting}
mdyn=closureDynamics(net,maxdeg,method,funname,symParameters)
\end{lstlisting}
\end{command}

\noindent
This \SDT{} function computes the exact (open) moment dynamics
\eqref{eq:exact} and then uses the moment closure technique specified
by the input parameter \lstinline{method} to compute approximate
(closed) moment dynamics of the form \eqref{eq:approx-1} or
\eqref{eq:approx-2}.
%
%\medskip 
%
As described in \cite{SinghHespanhaDec06}, the (exact) time-derivative
of the uncentered moment
\begin{align*}
  \mu^{(m)}&\eqdef\E[\psi(x_1,x_2,\dots,x_n)], &
  \psi(x_1,x_2,\dots,x_n)&\eqdef x_1^{m_1} x_2^{m_2}\cdots x_n^{m_n}
\end{align*}
can be obtained by computing
\begin{align}\label{eq:dynkin}
  \Deriv{\mu^{(m)}}{t} = \E[(L\psi)(x_1,x_2,\dots,x_k)]
\end{align}
where $(L\psi)(x_1,x_2,\dots,x_k)$ is an expression obtained by applying
to $\psi(\cdot)$ the generator of the Markov process that describes the
populations.
%
%\medksip
%
Regardless of the method used, the following is assumed in computing
the expected value in the left-hand side of \eqref{eq:dynkin}:
\begin{enumerate}
\item All replacement rules specified in the \lstinline{rules}
  section of the \texttt{.net} file are applied to
  $(L\psi)(x_1,x_2,\dots,x_k)$ \emph{before} taking the expected value
  in the right-hand side of \eqref{eq:dynkin}.
  
\item For every species $x_1$ that was declared
  \lstinline{deterministic}, it is assumed that
  \begin{align}\label{eq:deterministic}
    \E \big[ x_1 f(x_2,...x_k) \big] = \E[x_1] \E\big[ f(x_2,...x_k)
    \big]
  \end{align}
  for any function $f(.)$ of the remaining species.
  
\item For every species $x_1$ that was declared \lstinline{boolean},
  it is assumed that $x_1^n = x_1$, $\forall n\geq  1$. This introduces no
  error as long as the number of molecules of $x_1$ is indeed
  restricted to the set $\{0,1\}$.
\end{enumerate}

\paragraph{Inputs:}
\begin{itemize}
\item \lstinline{net} is a structure that describes the network of
  chemical reactions. It is typically obtained from a \texttt{.net}
  file using \lstinline{net=readNet(filename)}.
  
\item \lstinline{maxdeg} is an integer that specifies the largest degree
  for the uncentered moments in $\mu$. If \lstinline{maxdeg} is not an
  integer, then boolean variables are not taken into account for the degree
  of a moment. In this case, more moments are included in $\mu$ for the
  same value of \lstinline{maxdeg}.
  
\item \lstinline{method} (optional) is a string that specifies the
  method that should be used for moment closure.  This parameter is
  optional, and in its absence the first method below is used. The
  following methods are currently recognized:
  \begin{itemize}
  \item \lstinline{'dm'} or \lstinline{'derivativematching'} for
    moment closure obtained by matching the derivatives of the exact
    and approximate moment dynamics~\cite{SinghHespanhaDec06}. See
    Section~\ref{se:derivative-matching}.

  \item \lstinline{'zc'} or \lstinline{'zerocumulants'} for moment
    closure obtained by assuming that the high-order cumulants
    corresponding to all unknown moments are equal to zero. See
    Section~\ref{se:zero-cumulants}.

    For second order closure and elementary reactions, this corresponds
    to the technique described in \cite{UribeVerghese06}. However,
    this function can be used for closures of any order and for
    reactions with more than two reactants (non-elementary).

  \item \lstinline{'ld'} or \lstinline{'lowdispersion'} for moment
    closure obtained by assuming that the high-order normalized
    centered moments corresponding to all unknown moments are equal to
    zero. See Section~\ref{se:low-dispersion}.

    For second order closure and elementary reactions, this
    corresponds to the technique described in \cite{UribeVerghese06}.
    However, this function can be used for closures of any order and
    for reactions with more than two reactants (non-elementary).

  \item \lstinline{'qd'} or \lstinline{'quasideterministic'} for
    moment closure obtained by assuming that the high-order
    quasi-deterministic normalized centered moments corresponding to
    all unknown moments are equal to zero.

    For this approximation, the closed dynamics are of the form
    \eqref{eq:approx-2}, where $\phi$ is the solution to the
    determinstic dynamics (in molecule counts).  See
    Section~\ref{se:quasi-deterministic}.

  \item \lstinline!{'lna','Volume'}! or
    \lstinline!{'vankampen','Volume'}! or
    \lstinline!{'linearnoiseapproximation','Volume'}! for Van Kampen's
    linear noise approximation \cite[Chapter X]{Kampen01}. The second
    element of the cell specifies the variable to be used as the
    volume. This variable should have be declared in the
    \lstinline{parameters} section of the \texttt{.net} file.  To
    achieve moment closure, one considers the limit as this variable
    converges to infinity. For this limit to result in a closed set of
    moment equations, the all reactions should be elementary and their
    rates should depend on the volume as follows:
    \begin{align*}
      \emptyset  & \xrightarrow{\text{rate } = c\text{ \lstinline{Volume}}} * &
      X  & \xrightarrow{\text{rate } = c X} *  & 
      2 X  & \xrightarrow{\text{rate } = c\frac{X(X-1)}{\text{ \lstinline{Volume}}}} *  &
      X+Y  & \xrightarrow{\text{rate } = c\frac{X Y}{\text{ \lstinline{Volume}}}} * 
    \end{align*}
    For this approximation, the closed dynamics are of the form
    \eqref{eq:approx-2}, where $\phi$ is the solution to the
    determinstic dynamics (in concentrations). See
    Section~\ref{se:lna}.

  \item \lstinline{'zv'} or \lstinline{'zerovariance'} or
    \lstinline{'zv'} for moment closure assuming a negligible variance
    for the populations, i.e., assuming that the equality
    \eqref{eq:deterministic} holds for all species.
  \end{itemize}

\item \lstinline{funname} (optional) is a string containing the
  filename used to create an m-file that computes the left-hand side
  of the differential equations in \eqref{eq:approx-1} or
  \eqref{eq:approx-2}.  For approximate dynamics of the form
  \eqref{eq:approx-2}, the $\mu$ components of the state appear above
  the $\phi$ components. 
%
  When the function created by \lstinline{closureDynamics()} is called
  with no arguments, it outputs the names of the state variables (as a
  string array).
  
  The function created can be used as input to \lstinline{ode23s()} to
  simulate the moment dynamics (see Section~\ref{se:simulation} and
  the example in Section~\ref{se:quick}).
  
\item \lstinline{symParameters} (optional) is a vector of symbolic
  variables corresponding to parameters that appear in the rate
  expressions (e.g., rate constants).
  
  This vector is used only in the creation of the function
  \lstinline{funname} and it allows inclusion of additional inputs to
  this function that may be needed to specify numerical values for
  parameters that appear in the moment dynamics.  These values will
  override any default values specified in \lstinline{net.parameter}.
\end{itemize}


\paragraph{Outputs:}

\begin{itemize}
\item \lstinline{mdyn} is a structure characterizing the exact moment
  dynamics in \eqref{eq:exact} and the approximate moment dynamics in
  \eqref{eq:approx-1} or \eqref{eq:approx-2} (see
  Section~\ref{se:mdyn}). 
\end{itemize}

\paragraph{Under the hood:}
The function \lstinline{closureDynamics()} essentially combines the
functionality of the three lower-level functions:
\lstinline{momentDynamics()}, \lstinline{momentClosure()}, and
\lstinline{sym2mfile()}. Documentation for these functions is provided in
the corresponding m-files.

\section{The \lstinline{mdyn} structure}\label{se:mdyn}

The \lstinline{mdyn} structure characterizes the exact and approximate
moment dynamics, including
\begin{enumerate}
\item the entries of the $\mu$ vector in the exact and approximate
  moment dynamics (\lstinline{Mu} field),
\item the exact (open) moment dynamics in \eqref{eq:exact}
  (\lstinline{dotMu} field),
\item the approximate (closed) moment dynamics in \eqref{eq:approx-1}
  (\lstinline{approxDotMu} field) or in \eqref{eq:approx-2}
  (\lstinline{approxDotMu} and \lstinline{dotPhi} fields), and 
\item set of formatting rules to produce the moment dynamics in \LaTeX{}
  format (\lstinline{texrules} field).
\end{enumerate}
A detailed explanation of the key fields of this structure
follows. This structure may contain additional fields not described
here, which are used internally by \SDT{}.

\paragraph{Exact moment dynamics}
The following fields describe the exact moment dynamics in
\eqref{eq:exact}:
\begin{itemize}

\item \lstinline{Mu} is a structure describing the entries of the
  vector $\mu$ in \eqref{eq:exact} and the vector $\nu$ in
  \eqref{eq:approx-1} or \eqref{eq:approx-2-nu}. The following table
  describes the key fields of this structure for a network of chemical
  reactions with species \lstinline{X1,X2,X3,...}:
  \begin{center}\small
    \begin{tabular}{p{.6in}p{2in}p{3.1in}}\hline
      field & type & value of row corresponding to moment $E[ X1^{m1} X2^{m2} X3^{m3}\cdots]$\\\hline
      \lstinline!Mu.sym! & column vector of symbolic variables & 
      \lstinline!mu_X1m1_X2m2_X3m3...!\\
      \lstinline!Mu.mon! & column vector of symbolic expressions &
      \lstinline!X1^m1*X2^m2*X3^m3...!\\
      \lstinline!Mu.ndx!& matrix of integers with one column per
      species (in the order they were
      declared in the \texttt{.net} file) & 
      \lstinline!m1,m2,m3,...!\\
      \lstinline!Mu.x0! & column vector of doubles & initial value from the \texttt{.net} file\\\hline
    \end{tabular}
  \end{center}
    
\item \lstinline{barMu} is a structure describing the entries of the
  vector $\bar\mu$ in \eqref{eq:exact}. The following table describes
  the key fields of this structure for a network of chemical reactions
  with species \lstinline{X1,X2,X3,...}:
  \begin{center}\small
    \begin{tabular}{p{.6in}p{2in}p{3.1in}}\hline
      field & type & value of row corresponding\\&& to moment $E[ X1^{m1} X2^{m2} X3^{m3}\cdots]$\\\hline
      \lstinline!barMu.sym! & column vector of symbolic variables & 
      \lstinline!mu_X1m1_X2m2_X3m3...!\\
      \lstinline!barMu.mon! & column vector of symbolic expressions &
        \lstinline!X1^m1*X2^m2*X3^m3...!\\
      \lstinline!barMu.ndx!& matrix of integers with one column per
      species (in the order they were
      declared in the \texttt{.net} file) & 
      \lstinline!m1,m2,m3,...!\\\hline
    \end{tabular}
  \end{center}
  
\item \lstinline{dotMu} is a structure describing the exact
  moment dynamics \eqref{eq:exact}. It contains the following fields:
  \begin{center}\small
    \begin{tabular}{p{.6in}p{2in}p{3.1in}}\hline
      field & type & value\\\hline
      \lstinline!dotMu.A! & matrix of symbolic expressions & matrix $A$ in \eqref{eq:exact}\\
      \lstinline!dotMu.B! & matrix of symbolic expressions & matrix $B$ in \eqref{eq:exact}\\
      \lstinline!dotMu.sym! & column vector of symbolic expressions&whole right-hand side of
    \eqref{eq:exact}\\\hline
    \end{tabular}
  \end{center}
  The right-hand side of \eqref{eq:exact} can thus be obtained using
  any one of the following two symbolic expressions:
\begin{lstlisting}
   dotMu.A * Mu.sym + dotMu.B * barMu.sym
   dotMu.sym
\end{lstlisting}
\end{itemize}

\paragraph{Approximate moment dynamics}

The following fields describe the approximate moment dynamics in
\eqref{eq:approx-1} or \eqref{eq:approx-2-nu}, depending on the moment
closure technique used. It contains the following fields:
\begin{itemize}
\item \lstinline{approxDotMu} is a structure describing the right-hand
  sides of \eqref{eq:approx-1} or \eqref{eq:approx-2-nu}.  It contains
  the following fields:
  \begin{center}\small
    \begin{tabular}{p{1in}p{1.8in}p{2.9in}}\hline
      field & type & value\\\hline
      \lstinline!approxDotMu.barMu! & column vector of symbolic expressions & 
      approximate value of $\bar\mu$ in \eqref{eq:exact}\\
      \lstinline!approxDotMu.sym! & column vector of symbolic expressions & 
      right-hand sides of \eqref{eq:approx-1} or \eqref{eq:approx-2-nu}\\
      \lstinline!approxDotMu.A! & matrix of symbolic expressions &
      Jacobian of \lstinline!approxDotMu.sym! with respect to \lstinline!Mu.sym! \\
      \lstinline!approxDotMu.c! & vector of symbolic expressions &
      \lstinline!approxDotMu.sym-approxDotMu.A*Mu.sym!\\
      \hline
    \end{tabular}
  \end{center}
  The right-hand side of \eqref{eq:approx-1} or
  \eqref{eq:approx-2-nu} can thus be obtained using
  any one of the following three symbolic expressions:
\begin{lstlisting}
   dotMu.A * Mu.sym + dotMu.B * approxDotMu.barMu
   approxDotMu.sym
   approxDotMu.A * Mu.sym + approxDotMu.c
\end{lstlisting}
  The last representation is especially useful for the
  \lstinline!'qd'! and \lstinline!'lna'! moment closure methods
  because in this case the matrices \lstinline!approxDotMu.A! and
  \lstinline!approxDotMu.c! do not depend on the entries of
  \lstinline!Mu.sym!  [cf.~\eqref{eq:approx-3}].

  Currently, the moment closure method \lstinline!'lna'! does not return
  the field \lstinline!approxDotMu.barMu!

\end{itemize}
The following fields are only needed to describe the approximate
moment dynamics in \eqref{eq:approx-2}.
\begin{itemize}
\item \lstinline{Phi} is a structure describing the entries of the
  vector $\phi$ in \eqref{eq:approx-2-phi}. The following table describes
  the key fields of this structure for a network of chemical reactions
  with species \lstinline{X1,X2,X3,...}:
  \begin{center}\small
    \begin{tabular}{p{.6in}p{2in}p{3.1in}}\hline
      field & type & value of row corresponding to moment $E[ X1^{m1} X2^{m2} X3^{m3}\cdots]$\\\hline
      \lstinline!Phi.sym! & column vector of symbolic variables & 
      \lstinline!phi_X1m1_X2m2_X3m3...!\\
      \lstinline!Phi.mon! & column vector of symbolic expressions &
      \lstinline!X1^m1*X2^m2*X3^m3...! (\lstinline!'qd'! method)\\
      &&\lstinline!(X1/V)^m1*(X2/V)^m2*(X3/V)^m3...! (\lstinline!'lna'! method)\\
      \lstinline!Phi.x0! & column vector of doubles & initial value from the \texttt{.net} file\\\hline
    \end{tabular}
  \end{center}

\item \lstinline{dotPhi} is a vector of symbolic expressions with the
  derivative of $\phi$ in \eqref{eq:approx-2-phi}.
\end{itemize}
The following additional fields are only returned by the \lstinline{'lna'}
moment closure method:
\begin{itemize}
\item \lstinline{Xi} is a structure describing the entries of the
  first-order ``perturbation'' vector $\xi$ in \eqref{eq:xi-1}. The
  following table describes the key fields of this structure for a
  network of chemical reactions with species \lstinline{X1,X2,X3,...}:
  \begin{center}\small
    \begin{tabular}{p{.6in}p{2in}p{3.1in}}\hline
      field & type & value of row corresponding to\\&& $\E\big[ \big(\frac{X1-V \phi_1}{V^\frac{1}{2}}\big)^{m1}
      \big(\frac{X2-V\phi_2}{V^\frac{1}{2}}\big)^{m2} 
      \big(\frac{X3-V\phi_3}{V^\frac{1}{2}}\big)^{m3}\ldots\big]$\\\hline
      \lstinline!Xi.sym! & column vector of symbolic variables & 
      \lstinline!xi_X1m1_X2m2_X3m3...!\\
      \lstinline!Xi.mon! & column vector of symbolic expressions &
      \lstinline!((X1-V*phi_X1)/V^(1/2))^m1*...!\\
      \lstinline!Xi.mu2xi!& column vector of symbolic expressions& 
      expression for $\E\big[ \big(\frac{X1-V \phi_1}{V^\frac{1}{2}}\big)^{m1}
      \big(\frac{X2-V\phi_2}{V^\frac{1}{2}}\big)^{m2} 
      \big(\frac{X3-V\phi_3}{V^\frac{1}{2}}\big)^{m3}\ldots\big]$ as a
      function of the entries of $\mu$ and $\phi$\\
      \lstinline!Xi.xi2mu! & column vector of symbolic expressions& 
      expression for the uncentered moment $\E\big[
      X1^{m1}X2^{m2}X3^{m3}\ldots\big]$ in $\mu$ as a
      function of the entries of $\xi$ and $\phi$\\
      \lstinline!Xi.xi2barMu! & column vector of symbolic expressions& 
      expression for the uncentered moment $\E\big[
      X1^{m1}X2^{m2}X3^{m3}\ldots\big]$ in $\bar \mu$ as a
      function of the entries of $\xi$ and $\phi$\\\hline
    \end{tabular}
  \end{center}
    
  
\item \lstinline{dotXi} is a symbolic vector with the exact
  derivative of $\xi$ as in \eqref{eq:xi-1}.
  
\item \lstinline{approxDotXi} is a structure describing the
  approximate derivative of $\xi$ in the right-hand side of
  \eqref{eq:xi-2}. It contains the following fields:
  \begin{center}\small
    \begin{tabular}{p{1in}p{1.8in}p{2.9in}}\hline
      field & type & value\\\hline
      \lstinline!approxDotXi.sym! & column vector of symbolic variables & 
      right-hand side of \eqref{eq:xi-2}\\
      \lstinline!approxDotXi.A! & matrix of symbolic expressions &
      Jacobian of \lstinline!approxDotXi.sym! with respect to \lstinline!Xi.sym! \\
      \lstinline!approxDotXi.c! & vector of symbolic expressions &
      \lstinline!approxDotXi.sym-approxDotXi.A*Xi.sym!\\
      \hline
    \end{tabular}
  \end{center}
  The right-hand side of \eqref{eq:xi-2} can thus be obtained using
  any one of the following two symbolic expressions:
\begin{lstlisting}
   approxDotXi.sym
   approxDotXi.A * Xi.sym + approxDotXi.c
\end{lstlisting}
  The last representation is especially useful because the matrices
  \lstinline!approxDotXi.A! and \lstinline!approxDotXi.c! do not
  depend on the entries of \lstinline!Xi.sym!
  [cf.~\eqref{eq:approx-3}].

\end{itemize}  

\paragraph{\LaTeX{} formatting} The following field is used to produced
\LaTeX-formatted versions of any symbolic expression involving
reaction parameters, moments, or exact and approximate moment
dynamics.  These rules should be used by the \SDT{} function
\lstinline{mylatex()} described in Section~\ref{se:utilities}.
\begin{itemize}
\item \lstinline{texrules} is a cell array of strings with one
  row per transformation rule and two columns:
  \begin{itemize}
  \item the first column contains a regular expression,
    following the syntax of \matlab{}'s \lstinline{regexp()}
  \item the second column contains the string to replace the
    regular expression.
  \end{itemize}
\end{itemize}

\chapter{Monte Carlo Simulations}

One can run Monte Carlo simulations of a network of chemical reactions
specified by a \texttt{.net} file within \matlab{} using \SDT or
in C++ using \textsc{StochKit}. 

\medskip
\noindent
\textbf{Attention!} For both simulation methods, 

\begin{enumerate}
\item the \lstinline{post_reaction_counts} in the \texttt{.net} file
  must correspond to increments/decrements; and
\item molecule counts are assume constant between reactions and
  therefore the \lstinline{derivatives} section of the \texttt{.net}
  file is ignored.
\end{enumerate}


\section{Using \SDT}


\begin{command}{quadPropensities()}
\begin{lstlisting}
[Q,b,c,s,x0]=quadPropensities(net)
\end{lstlisting}
\end{command}

\noindent
The function
\lstinline{quadPropensities()} describes a network of chemical reactions with quadratic propensity functions. The outputs of this function are used by \lstinline{sampledSSA()}
to run Monte Carlo simulations of the network of chemical reactions.

\paragraph{Inputs:}
\begin{itemize}
\item \lstinline{net} is a structure that describes the network of
  chemical reactions. It is typically obtained from a \texttt{.net}
  file using \lstinline{net=readNet(filename)}.
\end{itemize}

\paragraph{Outputs:}
\begin{itemize}
\item \lstinline{Q} is a matrix with the quadratic terms for the
  propensity functions. Each reaction corresponds to a (\# species) by
  (\# species) square matrix. These matrices are stacked on top of
  each other so the size of \lstinline{Q} is (\# reactions $\times$ \#
  species) by (\# species). 
  
\item \lstinline{b} is a matrix with the linear terms for the
  propensity functions.  Each reaction corresponds to a vector of
  length (\# species). These vectors are stacked on top of each other
  so size of \lstinline{b} is (\# reactions) by (\# species).
  
\item \lstinline{c} is a vector with the constant terms for the
  propensity functions.  The constants for each reaction are stacked
  on top of each other so size of \lstinline{c} is (\# reactions) by 1.
  
\item \lstinline{s} is a matrix with the stoichiometry for the
  reactions.  Each reaction corresponds to a vector of
  length (\# species). These vector are stacked
  side by side so size of \lstinline{s} is (\# species) by (\# reactions).
  
\item \lstinline{x0} is a vector with the initial conditions in the
  \texttt{.net}  file. The size of \lstinline{x0} is  (\# species) by
  1.
\end{itemize}

\begin{command}{sampledSSA()}
\begin{lstlisting}
[X,Xmean,Xstd,XmeanCI,XstdCI] = sampledSSA(Q,b,c,s,x0,nMC,Ts,CI)
\end{lstlisting}
\end{command}

\noindent
The function
\lstinline{sampledSSA} runs multiple Gillespie's stochastic Simulation Algorithms (SSA) for the network of chemical reactions with propsenties specified by the inputs \lstinline{Q},
\lstinline{b},
\lstinline{c} and stoichiometry specified by \lstinline{s}.

This function returns all the sample paths at a set of sample times,
as well as estimates for the means and standard deviations at the
sample times.

\paragraph{Inputs:}
\begin{itemize}
\item \lstinline{Q} is a matrix with the quadratic terms for the
  propensity functions. Each reaction corresponds to a (\# species) by
  (\# species) square matrix. These matrices are stacked on top of
  each other so the size of \lstinline{Q} is (\# reactions $\times$ \#
  species) by (\# species). 
  
\item \lstinline{b} is a matrix with the linear terms for the
  propensity functions.  Each reaction corresponds to a vector of
  length (\# species). These vectors are stacked on top of each other
  so size of \lstinline{b} is (\# reactions) by (\# species).
  
\item \lstinline{c} is a vector with the constant terms for the
  propensity functions.  The constants for each reaction are stacked
  on top of each other so size of \lstinline{c} is (\# reactions) by 1.
  
\item \lstinline{s} is a matrix with the stoichiometry for the
  reactions.  Each reaction corresponds to a vector of
  length (\# species). These vector are stacked
  side by side so size of \lstinline{s} is (\# species) by (\# reactions).
  
\item \lstinline{x0} is a vector with the initial conditions in the
  \texttt{.net}  file. The size of \lstinline{x0} is  (\# species) by
  1.

\item \lstinline{nMC} is the number if Monte Carlo simulations to run.

\item \lstinline{Ts} is a vector with the desired sample times for the
  output.

  The sampled paths are computed exactly at all times, but their
  values are only returned at the time in the vector \lstinline{Ts}.

\item \lstinline{CI} is an optional input with the desired percentage
  for the confidence intervals for the mean and standard deviation. If
  not specified, 95\% is used.
\end{itemize}

\paragraph{Outputs:}
\begin{itemize}
\item \lstinline{X} is matrix with the molecule counts at the sample
  times. The size of \lstinline{X} is (\# species) by (\# simulations)
  by (\# of sample times).

\item \lstinline{Xmean} is a matrix with the mean molecule counts at
  the sample times. The size of \lstinline{Xmean} is (\# species) by
  (\# of sample times).

\item \lstinline{Xstd} is a matrix with the standard deviation
  molecule counts at the sample times.  The size of \lstinline{Xstd}
  is (\# species) by (\# of sample times)

\item \lstinline{XmeanCI} is a matrix with the confidence interval for
  the \lstinline{Xmean}. The size of \lstinline{XmeanCI} is (\#
  species) by (\# of sample times) by 2.

  The computation of \lstinline{XmeanCI} assumes that the central
  limit theorem is valid to determine the distribution of the mean

\item \lstinline{XstdCI} is a matrix with the confidence interval for
  the Xstd.  The size of XstdCI is (\# species) by (\# of sample
  times) by 2.

  The computation of \lstinline{XstdCI} assumes a normal distribution.
\end{itemize}



\section{Using  \textsc{StochKit}}

\begin{command}{net2stochKit()}
\begin{lstlisting}
net2stochKit(net,filename,x0,symParameters,valueParameters)
\end{lstlisting}
\end{command}

\noindent
The function \lstinline{net2stochKit()} automatically creates a C++ file that
can be used to create a \textsc{StochKit} executable to run Monte Carlo
simulations of a network of chemical reactions.

\medskip

\textbf{Attention!} All the \lstinline{post_reaction_counts} in the
\texttt{.net} file must correspond to increments/decrements. If this
is not the case the C++ code will not compile.

\paragraph{Inputs and Outputs:}
\begin{itemize}
\item \lstinline{net} is a structure that describes the network of
  chemical reactions. It is typically obtained from a \texttt{.net}
  file using \lstinline{net=readNet(filename)}.
  
\item \lstinline{filename} is a string containing the filename of the
  C++ file to be created (without the \lstinline{.cpp} extension).
  
\item \lstinline{x0} is a vector of initial populations for the Monte
  Carlo simulations with as many entries as the number of chemical
  species.
  
\item \lstinline{symParameters} (optional) is a vector of symbolic
  parameters that appear in the rate expressions (e.g., rate
  constants).
  
\item \lstinline{valueParameters} (optional) is a vector of numerical
  values for the symbolic parameters in symParameters.
\end{itemize}
The function \lstinline{net2stochKit()} returns no output, but creates a
C++ \textsc{StochKit} ``ProgramDefinition'' file. To learn how to use
this file, please consult the StochKit user guide \cite{StochKit}. An
example Makefile and two ``main'' C++ files are provided as examples.
However, these will only work if (i) \textsc{StochKit} has been
successfully installed and (ii) the variable \lstinline{CSE_CPP_HOME} in
the Makefile points to the \textsc{StochKit} directory.

\medskip

An error is returned if equations for continuous evolution are
specified in the \lstinline{derivatives} section because
\textsc{StockKit} cannot simulate reaction networks with continuous
variations in the populations.

\chapter{Simulation and plots}\label{se:simulation}

The m-file produced by \lstinline{closureDynamics()} can be passed to
an ODE solver such as \lstinline{ode23s()} to solve the approximate
moment dynamics. The solution to the ODE can then be used to compute
and plot centered moments using \lstinline{getCMoments()} and
\lstinline{plotCMoments()}, respectively. It can also be used to
obtain an approximate probability distribution using
\lstinline{getDistribution()}.


\medskip

\begin{command}{ode23s()}
\begin{lstlisting}
[t,mu]=ode23s(@(t,x)funname(x,parameter_list),[0,Tmax],x0);
\end{lstlisting}
\end{command}

\noindent
This (standard) \matlab{} function can be used to solve the approximate
moment dynamics.

\paragraph{Inputs:}
\begin{itemize}
\item \lstinline{funname} is a string containing the filename of the
  m-file created by \lstinline{closureDynamics()} to compute the
  left-hand side of \eqref{eq:approx-1} or \eqref{eq:approx-2}
  (without the \lstinline{.m} extension).
  
\item \lstinline{parameter_list} is a list of parameter values to be
  passed to \lstinline{funname()}, as specified by the input parameter
  \lstinline{symParameters} to the function \lstinline{closureDynamics()}.
  
\item \lstinline{Tmax} is the time at which the simulation should
  terminate.
  
\item \lstinline{x0} is the initial condition for the moment dynamics.

  For moment dynamics of the form \eqref{eq:approx-1}, \lstinline{x0}
  must have the size of $\mu$ and to obtain the initial the conditions
  in the \texttt{.net} file one would choose
  \lstinline{x0=mdyn.Mu.x0}, with \lstinline{mdyn} returned by
  \lstinline{closureDynamics()}.

  For moment dynamics of the form \eqref{eq:approx-2}, \lstinline{x0}
  must have the size of $\matt{\mu' & \phi'}'$ and to obtain the initial the
  conditions in the \texttt{.net} file one would choose
  \lstinline{x0=[mdyn.Mu.x0;mdyn.Phi.x0]}, with \lstinline{mdyn}
  returned by \lstinline{closureDynamics()}.
  
  \textbf{Attention!} Starting with initial conditions for which some
  species have exactly zero molecules can be problematic when using
  moment closure based on derivative matching (because this may lead
  to divisions by zero or errors in the ODE solver). To avoid this,
  one can initialize the system with small but positive population
  (say .1 molecules).
\end{itemize}

The function \lstinline{ode23s()} accepts several other parameters (see
\lstinline{help ode23s}). Other ODE solvers can also be used, however a
stiff solver is generally preferable.

\paragraph{Outputs:}

\begin{itemize}
\item \lstinline{t} is a column vector with the time instants at which the
  solution was computed.
  
\item \lstinline{mu} is a time series of $\mu$ for moment dynamics of
  the form \eqref{eq:approx-1} or a time series of $\matt{\mu' & \phi'}'$
  for moment dynamics of the form \eqref{eq:approx-2}.  This vector
  can be plotted using the following command
  \begin{center}\codesize
    \begin{tabular}{|c|}\hline
      \begin{minipage}[c]{5.25in}
\begin{lstlisting}
plot(t,mu)
legend(funname(),'location','best','interpreter','none')
\end{lstlisting}
      \end{minipage}\\\hline
    \end{tabular}
  \end{center}  
  where \lstinline{funname()} is the m-file created by
  \lstinline{closureDynamics()}. Recall that when the function created
  by \lstinline{closureDynamics()} is called without arguments it
  returns the names of the variables in $\mu$, which can be used to
  label the plot.  The option \lstinline{'interpreter','none'} to the
  command \lstinline{legend()} prevents \matlab{} from interpreting
  underscores as subscripts.
\end{itemize}

\begin{command}{getCMoments()}
\begin{lstlisting}
[average,stddev]=getCMoments(net,mdyn,mu,symExpression,assumezero)
\end{lstlisting}
\end{command}

\noindent
Given a time series $\mu(t)$ of uncentered moments, typically computed
using \lstinline{ode23s()}, this \SDT{} function computes the means
and standard deviations of a given vector of expressions involving the
populations of the different species.

\paragraph{Inputs:}
\begin{itemize}
\item \lstinline{net} is a structure that describes the network of
  chemical reactions. It is typically obtained from a \texttt{.net}
  file using \lstinline{net=readNet(filename)}.

\item \lstinline{mdyn} is a structure characterizing the exact moment
  dynamics in \eqref{eq:exact} and/or the approximate moment dynamics
  in \eqref{eq:approx-1} or \eqref{eq:approx-2}, typically computed
  using \lstinline{closureDynamics()}.

\item \lstinline{mu} is a matrix containing a time series of $\mu$ or
  $\matt{\mu' & \phi'}'$, typically computed using
  \lstinline{ode23s()}.  Each row of \lstinline{mu} corresponds to a
  different time instant and each column to a different moment.
  
\item \lstinline{symExpression} is an array of symbolic expressions
  whose means and standard deviations will be computed. The function
  returns \lstinline{NaN} if the computation of a mean and/or
  standard deviation cannot be computed exactly using the moments
  available (except when \lstinline{assumezero} is nonzero, see below).
  
  The replacement rules in \lstinline{net.substitutionrule} will be
  applied to the symbolic expression (and to its square in order to
  compute the variance) before taking expectations. The boolean nature
  of the variables declared as such, will be taken into account.
  
  \textbf{Attention!} In the current implementation of this function,
  the deterministic nature of the variables declared as such will be
  ignore. This means, e.g., that \lstinline{NaN} will be generated if
  one asks for the expected value of \texttt{X\^{}2} and
  \lstinline{mu} only contains the expected value of \lstinline{X},
  \emph{even if \lstinline{X} was declared deterministic}.

\item \lstinline{assumezero} is an optional boolean variable that,
  when nonzero, instructs the function to assume zero all moments that
  are not available to compute the mean and/or standard deviations.
\end{itemize}

\paragraph{Outputs:}
\begin{itemize}
\item \lstinline{average} is a time series of the expected value of the
  vector \lstinline{symExpression}, with one time instant per row and one
  element of the \lstinline{symExpression} per column.
  
\item \lstinline{stddev} (optional) is a time series of the standard
  deviation of the vector \lstinline{symExpression}, with one time
  instant per row and one element of the \lstinline{symExpression} per
  column.
  
  \textbf{Attention!} For populations with low stochasticity (i.e.,
  for which standard deviations are much smaller than the means), the
  moment closure approximation errors can lead to negative variances.
  Such populations should be declared as \lstinline{deterministic} in the
  \texttt{.net} file.
\end{itemize}


\begin{command}{plotCMoments()}
\begin{lstlisting}
[handles,average,stddev]=plotCMoments(net,mdyn,t,mu,symExpression_1,style_1,symExpression_2,style_2,...)
\end{lstlisting}
\end{command}

\noindent
Given a time series $\mu(t)$ of uncentered moments, typically computed
using \lstinline{ode23s()}, this \SDT{} function computes and plots
the means and standard deviations of a given vector of expressions
involving the populations of the different species. Internally, this
function uses \lstinline{getCMoments()} to compute means and standard
deviations.

\paragraph{Inputs:}
\begin{itemize}
\item \lstinline{net} is a structure that describes the network of
  chemical reactions. It is typically obtained from a \texttt{.net}
  file using \lstinline{net=readNet(filename)}.

\item \lstinline{mdyn} is a structure characterizing the exact moment
  dynamics in \eqref{eq:exact} and/or the approximate moment dynamics
  in \eqref{eq:approx-1} or \eqref{eq:approx-2}, typically computed
  using \lstinline{closureDynamics()}.

\item \lstinline{mu} is a matrix containing a time series of $\mu$ or
  $\matt{\mu' & \phi'}'$, typically computed using
  \lstinline{ode23s()}.  Each row of \lstinline{mu} corresponds to a
  different time instant and each column to the different moment.
  
\item \lstinline{symExpression_1,symExpression_2,...} are symbolic
  expressions whose means and standard deviations will be plotted. The
  means of each $\lstinline{symExpression_i}$ will be plotted with the
  style defined by \lstinline{style_i} and dotted lines will be used
  to depict the mean plus/minus one standard deviation.
  
  This function uses \lstinline{getCMoments()} so the reader is referred
  to the documentation of that function for details on the computation
  of means and standard deviations.
  
\item \lstinline{style_1,style_2,...} are character strings defining
  the style of the line to be used, as in \matlab{}'s \lstinline{plot()}
  command.
\end{itemize}

\paragraph{Outputs:}
\begin{itemize}
\item \lstinline{handles} is a vector with the handles of all the
  lines plotted.  These handles are useful, for example to produce
  legends, as in
\begin{center}\codesize
  \begin{tabular}{|c|}\hline
    \begin{minipage}[c]{5.25in}
\begin{lstlisting}
h=plotCMoments(net,monMu,t,mu,'E','b-','St','g-','Et','k-');
legend(h(1:3:end),'E[E]\pm{}Std[E]','E[St]\pm{}Std[St]','E[Et]','location','best');
\end{lstlisting}
    \end{minipage}\\\hline
  \end{tabular}
\end{center}  
where by using only one out of each three handles, one skips legends
for the lines corresponding to the standard deviations.

\item \lstinline{average} is a time series with the expected values of
  all expressions, with one time instant per row and one expression
  per column.
  
\item \lstinline{stddev} is a time series with the standard deviations
  of all expressions, with one time instant per row and one expression
  per column.
\end{itemize}

\begin{command}{getDistribution()}
\begin{lstlisting}
pdf=getDistribution(net,mdyn,mu,symExpression,x,distribution,symBoolean)
\end{lstlisting}
\end{command}

\noindent
Given a vector $\mu$ of uncentered moments, this \SDT{} function
estimates the distribution of a given expression involving the
populations of the different species. Internally, this function uses
\lstinline{getCMoments()} to compute means and standard deviations.

\paragraph{Inputs:}
\begin{itemize}
\item \lstinline{net} is a structure that describes the network of
  chemical reactions. It is typically obtained from a \texttt{.net}
  file using \lstinline{net=readNet(filename)}.

\lstinline{mdyn} is a structure characterizing the exact moment
  dynamics in \eqref{eq:exact} and the approximate moment dynamics in
  \eqref{eq:approx-1} or \eqref{eq:approx-2}

\item \lstinline{mu} is a vector containing the value of $\mu$ or
  $\matt{\mu' & \phi'}'$ at a given time instant, typically computed
  using \lstinline{ode23s()}.

\item \lstinline{symExpression} is a symbolic expression whose
  distribution will be estimated.
  
  This function uses \lstinline{getCMoments()} to compute the mean and
  standard deviation from which the distribution is estimated. The
  reader is referred to the documentation of \lstinline{getCMoments()}
  for details on the computation of means and standard deviations.

\item \lstinline{x} is a vector containing the points at which the
  distribution will be computed.
  
\item \lstinline{distribution} is a character string specifying which
  type of distribution should be assumed. This parameter can take the
  following values:
  \begin{itemize}
  \item \lstinline{'normal'} when a normal distribution should be
    assumed.

  \item \lstinline{'lognormal'} when a lognormal distribution should be
    assumed.
    
  \item \lstinline{'binomial'} when a binomial distribution should be
    assumed.
    
  \item \lstinline{'mix_normal'} when it should be assumed that the
    distribution is normal, when conditioned to any one of the two
    values for a given boolean-valued expression \lstinline{symBoolean}.
    This will result in a convex combination (mixture) of two
    normal distributions.

  \item \lstinline{'mix_lognormal'} when it should be assumed that
    the distribution is lognormal, when conditioned to any one of the
    two values for a given boolean-valued expression
    \lstinline{symBoolean}. This will result in a convex combination
    (mixture) of two lognormal distributions.
  \end{itemize}
  
\item \lstinline{symBoolean} is a boolean-valued symbolic expression
  needed for any of the ``mixture'' distributions discussed above.
\end{itemize}  

\paragraph{Output:}
\begin{itemize}
\item \lstinline{pdf} is a vector containing the values of the
  probability density function (pdf) at the points in \lstinline{x}.
  \end{itemize}


\begin{command}{plotCompare()}
\begin{lstlisting}
[exact,approx,MC]=plotCompare(net,yPlots,xPlots,Tmax,varargin)
\end{lstlisting}
\end{command}

\noindent
This function produces an array of plots that can be used to compare
several moment closure techniques and Monte Carlo Simulations. 

\medskip

\emph{This function is currently undocumented but we provide an
  example of its use below. More examples can be found in the
  \lstinline{/examples} folder}.

\begin{lstlisting}
figure(1);clf
[exact,approx,MC]=plotCompare(net,2,5,Tmax,...
    {{'derivative matching',X1minmax},...                  % subfigure (2,5,1)
           {'MCsample',1,'X1','c-'},{'dm',2,'X1','b-'},{'dm',3,'X1','r-'}},...
    {{'zero cumulants',X1minmax},...                       % subfigure (2,5,2)
           {'zc',1,'X1','k-'},{'zc',2,'X1','b-'},{'zc',3,'X1','r-'}},...
    {{'low dispersion',X1minmax},...                       % subfigure (2,5,3)
           {'ld',1,'X1','k-'},{'ld',2,'X1','b-'},{'ld',3,'X1','r-'}},...
    {{'quasi deterministic',X1minmax},...                  % subfigure (2,5,4)
           {'qd',1,'X1','k-'},{'qd',2,'X1','b-'},{'qd',3,'X1','r-'}},...
    {{'distribution (Monte Carlo)',X1minmax},...           % subfigure (2,5,5)
           'MCdistribution',nMC,'X1'},...   
    {{'derivative matching',X2minmax},..                   % subfigure (2,5,6).
           {'MCsample',1,'X2','c-'},{'dm',2,'X2','b-'},{'dm',3,'X2','r-'}},...
    {{'zero cumulants',X2minmax},...                       % subfigure (2,5,7)
           {'zc',1,'X2','k-'},{'zc',2,'X2','b-'},{'zc',3,'X2','r-'}},...
    {{'low dispersion',X2minmax},...                       % subfigure (2,5,8)
           {'ld',1,'X2','k-'},{'ld',2,'X2','b-'},{'ld',3,'X2','r-'}},...
    {{'quasi deterministic',X2minmax},...                  % subfigure (2,5,9)
           {'qd',1,'X2','k-'},{'qd',2,'X2','b-'},{'qd',3,'X2','r-'}},...
    {{'distribution (Monte Carlo)',X2minmax},...           % subfigure (2,5,10)
           'MCdistribution',nMC,'X2'});
\end{lstlisting}
In this example,
\begin{enumerate}
\item \lstinline{net} is a structure that describes the network of
  chemical reactions. It is typically obtained from a \texttt{.net}
  file using \lstinline{net=readNet(filename)}.

\item \lstinline{2,5} specifies that an array of 2 by 5 subfigure
  should be draw. Each subfigure can have several plots. All subplots
  for the same figure are passed as a cell array.

\item \lstinline{Tmax} specifies the that all simulations should be
  produced in the interval \lstinline{[0,Tmax]} and that Monte
  Carlos distributions should be computed for the time
  \lstinline{Tmax}.

\item \lstinline{'derivative matching',X1minmax} specifies that the tile of
  the first subfigure should be \lstinline{'derivative matching'} and
  \lstinline{X1minmax} is a 2-vector that specifies the limits of
    the y-axis for this subfigure.

\item \lstinline{'MCsample',1,'X1','c-'} specifies a plot within a
  subfigure, containing 
  \lstinline{1} Monte Carlo sample path of the species \lstinline{X1},
  using the line style \lstinline{c-}.

\item \lstinline{'dm',2,'X1','b-'} specifies a plot within a
  subfigure, displaying
  the moment closure approximate dynamics using derivative matching
  with moments up to order \lstinline{2}. The plot with contain the
  mean of \lstinline{X1} using the line style \lstinline{b-}, and
  dotted lines specifying $\pm1$ standard deviation.

\item \lstinline{'MCdistribution',nMC,'X2'} specifies a plot within a
  subfigure, displaying an histogram
  with the distribution of \lstinline{X2} at time \lstinline{Tmax}
  obtained from  \lstinline{nMC} Monte Carlo simulations.
\end{enumerate}

\chapter{Utilities}\label{se:utilities}

\begin{command}{subsParameters()}
\begin{lstlisting}
new_expression=subsParameters(net,old_expression)
\end{lstlisting}
\end{command}

The function \lstinline{subsParameters()} takes a given expression and
replaces the parameters from a network of chemical reactions by their
default values.

\paragraph{Inputs and Outputs:}
\begin{itemize}
\item \lstinline{net} is a structure that describes the network of
  chemical reactions. It is typically obtained from a \texttt{.net}
  file using \lstinline{net=readNet(filename)}.
  
\item \lstinline{old_expression} is a vector of symbolic expressions
  where the parameters should be eliminated.
  
\item \lstinline{new_expression} is a vector of symbolic expressions
  obtained from \lstinline{old_expression} by replacing any parameters
  by their default values.
\end{itemize}

\begin{command}{mylatex()}
\begin{lstlisting}
str=mylatex(sym,texrules)
\end{lstlisting}
\end{command}

The function \lstinline{mylatex()} takes a symbolic expression and
returns it into a fairly readable \LaTeX{} syntax. This function allows
the latex expression to be transformed by a set of replacement rules,
which can be used to include the \LaTeX{} names specified in the
\texttt{.net} file.

\paragraph{Inputs and Outputs:}

When calling \lstinline{mylatex()} without inputs, this function returns
a set of macros needed to interpret the \LaTeX{} output. Otherwise\ldots
  
\begin{itemize}
\item \lstinline{sym} is a symbolic expression to be transformed into
  \LaTeX{} format.

\item \lstinline{texrules} is a cell array of strings with one row per
         transformation rule and two columns:
         \begin{itemize}
         \item the first column contains a regular expression,
           following the syntax of \matlab{}'s \lstinline{regexp()}
         \item the second column contains the string to replace the
           regular expression.
         \end{itemize}
         The transformation rules are applied sequentially and only
         once. Replacements are not applied to strings already replaced.

         \lstinline{texrules} are typically obtained by the \SDT{} command
         \lstinline{closureDynamcis()} or \lstinline{momentDynamics()}.

\item \lstinline{str} is a string in the \LaTeX{} output.
\end{itemize}

\begin{command}{momentEquilibrium()}
\begin{lstlisting}
[muEq,phiEq,xiEq]=momentEquilibrium(net,mdyn,option)
\end{lstlisting}
\end{command}

This \SDT{} function computes the equilibrium points of approximate
moment dynamics of the form \eqref{eq:approx-1} or \eqref{eq:approx-2}.

\paragraph{Inputs:}
\begin{itemize}
\item \lstinline{net} is a structure that describes the network of
  chemical reactions. It is typically obtained from a \texttt{.net}
  file using \lstinline{net=readNet(filename)}.

\item \lstinline{mdyn} is a structure describing the moment dynamics in
  \eqref{eq:approx-1} or \eqref{eq:approx-2}.  This structure is
  typically produced by \lstinline{closureDynamics()}.
    
\item \lstinline{numerical} is an optional boolean variable indicating
  whether or not the parameters should be replaced by the numerical
  values in \lstinline{net.parameter} prior to finding the equilibrium
  points.

\item \lstinline{option} is an optional string specifying options to be
  used in the computation of the equilibrium points. The following
  options are currently recognized:
  \begin{itemize}
  \item \lstinline{'none'} --- no option specified.

  \item \lstinline{'lna'} --- use Van Kampen's linear noise approximation
    (see~\cite[Chapter X]{Kampen01}). Assumes that the \lstinline{'lna'}
    moment closure method was used to produce \lstinline{mdyn}.

  \item \lstinline{'smallCV'} --- assume that the distribution of the
    populations is tightly concentrated around a mean value.
    Currently does not work for Van Kampen's linear noise
    approximation.
    \end{itemize}
\end{itemize}


\paragraph{Outputs:}

\begin{itemize}
\item \lstinline{muEq} is a cell array of symbolic vectors with the
  equilibrium points of the vector $\mu$ in \eqref{eq:approx-1} or
  \eqref{eq:approx-2}. Each entry of the cell array corresponds to
  one equilibrium point.
\end{itemize}
The following additional outputs are available with the options
\lstinline{'lna'} and \lstinline{'smallCV'} 
\begin{itemize}
\item \lstinline{phiEq} is a cell array of symbolic vectors with the
  equilibrium points of the deterministic vector $\phi$ around which
  the populations are concentrated.  Each entry of the cell array
  corresponds to one equilibrium point.

\item \lstinline{xiEq} is a cell array of symbolic vectors with the
  equilibrium points of the perturbations $\xi$ around $\phi$.  Each
  entry of the cell array corresponds to one equilibrium point.
\end{itemize}


\appendix
\chapter{Appendix}

\section{Installation}

The \lstinline{.tar} file posted contains the following three folders:
\begin{itemize}
\item \lstinline{stochdyntools/doc} --- Packet documentation (including this
  file).
\item \lstinline{stochdyntools/lib} --- The \SDT{} \matlab{} scripts described in
  this document.
  
\item \lstinline{stochdyntools/examples} --- Examples on the use of
  this set of tools (including the examples in
  Section~\ref{se:quick}).
\end{itemize}
To start using \SDT{} one needs to
\begin{enumerate}
\item ``Untar'' the \lstinline{.tar} file. In a PC this can be done using
  any of the usual ``unzip'' utilities and in Linux this can be done
  by executing the following command at the shell prompt:
  \begin{center}\codesize
    \begin{tabular}{|c|}\hline
      \begin{minipage}[c]{5.25in}
\begin{lstlisting}
tar xvf sotchdyntools.tar
\end{lstlisting}
      \end{minipage}\\\hline
    \end{tabular}
  \end{center}  

\item Add the \lstinline{stochdyntools/lib} folder to the \matlab{}
  path. This can be done by executing the following command at the
  \matlab{} prompt:
  \begin{center}\codesize
    \begin{tabular}{|c|}\hline
      \begin{minipage}[c]{5.25in}
\begin{lstlisting}
path('{BASE}stochdyntools/lib',path)
\end{lstlisting}
      \end{minipage}\\\hline
    \end{tabular}
  \end{center}  
where \texttt{\{BASE\}} should be replaced by the folder where the
\lstinline{.tar} file was ``untared.'' Most likely, you will want to add
the above command to your \lstinline{startup.m} \matlab{} script.
\end{enumerate}

If one wants to use the function \lstinline{net2stochKit()} to create
Monte Carlos simulations, \textsc{StochKit} must be installed
\cite{StochKit}. However, this is not needed to compute moment
dynamics.

\section{Algebraic derivations}

\section{Cumulants}\label{se:cumulants}

The multi-variable cumulants for the distribution of a random vector
$x$ are defined by the Taylor series of the cumulant-generating
function:
\begin{align*}
  g(\lambda)\eqdef \log\big(\E[e^{j \lambda \cdot x}]\big),
\end{align*}
where $\lambda\in\R^n$ and $\cdot$ denotes inner product. In particular, the
cumulant $\kappa^{(m)}$ associated with the vector of integers
$m\eqdef(m_1,m_2,\dots,m_n)$ is given by
\begin{align*}
  \kappa^{(m)}\eqdef 
  \PDeriv{^{m_1+m_2+\cdots+m_n} g(0)}{\lambda_1^{m_1}\lambda_2^{m_2}\cdots\lambda_n^{m_n}},
\end{align*}
which is one of the coefficients in the Taylor series of $g(\lambda)$ around
$\lambda=0$. The cumulant $\kappa^{(m)}$ can be computed by first expanding 
\begin{align*}
  \E[e^{j \lambda \cdot x}] 
  = \sum_{\ell=0}^{\sum m_i} \frac{1}{\ell!}\E\Big[(j\sum_i \lambda_i x_i)^\ell\Big]
  +\E[R(\lambda,x)],
\end{align*}
which results in 
\begin{align*}
  \E[e^{j\lambda \cdot x}] =\sum_{\sum \bar m_i\leq \sum m_i} c_{\bar m}(\lambda) \mu^{(\bar m)}
  +\E\Big[O\big((\lambda\cdot x)^{1+\sum m_i}\big)\big],
\end{align*}
where the summation is over moments $\mu^{(\bar m)}$ of order up to $\sum m_i$,
the $c_{\bar m}(\lambda)$ are appropriately selected constants that depend on $\lambda$,
and all terms in the remainder $R(\lambda,x)$ have powers in $\lambda$ or order
higher than $\sum m_i$. The cumulant $\kappa^{(m)}$ is then obtained from 
\begin{align*}
  \kappa^{(m)}
  %\eqdef\PDeriv{^{m_1+m_2+\cdots+m_n}g(0)}{\lambda_1^{m_1}\lambda_2^{m_2}\cdots\lambda_n^{m_n}}
  = \sum_{\sum \bar m_i\leq \sum m_i}  \mu^{(\bar m)}\PDeriv{^{m_1+m_2+\cdots+m_n}c_{\bar m}(0)}{\lambda_1^{m_1}\lambda_2^{m_2}\cdots\lambda_n^{m_n}}.
\end{align*}
Note that since the derivative is computed at $\lambda=0$, the terms in
$R(\lambda,x)$ will note appear because they only have powers in $\lambda$ or
order higher than $\sum m_i$. In fact, it is straightforward to show that
\begin{align*}
  \kappa^{(m)}
  = \mu^{(m)}+\sum_{\sum \bar m_i<\sum m_i}  \mu^{(\bar m)}\PDeriv{^{m_1+m_2+\cdots+m_n}c_{\bar m}(0)}{\lambda_1^{m_1}\lambda_2^{m_2}\cdots\lambda_n^{m_n}},
\end{align*}
which means that the cumulant $\kappa^{(m)}$ depends only on $\mu^{(m)}$ and
lower-order moments, so by setting $\kappa^{(m)}=0$ one obtains an
expression for $\mu^{(m)}$ as a function of lower-order moments. 


\section{Uncentered Versus Normalized Centered moments}\label{se:centered}

The uncentered moment $\mu^{(m)}$ can be expressed in terms of the
normalized centered moment $\eta^{(\bar m)}$ as follows
\begin{align*}
  \mu^{(m)} 
  &= \E[ x_1^{m_1} x_2^{m_2}\cdots x_n^{m_n} ] \\
  &= \E[x_1]^{m_1}\E[x_2]^{m_2}\cdots\E[x_n]^{m_n}
  \E\Big[\Big(1+\frac{x_1-\E[x_1]}{\E[x_1]}\Big)^{m_1}
  \Big(1+\frac{x_2-\E[x_2]}{\E[x_2]}\Big)^{m_2}\cdots
  \Big(1+\frac{x_n-\E[x_n]}{\E[x_n]}\Big)^{m_n}\Big]\\
  &=\E[x_1]^{m_1}\E[x_2]^{m_2}\cdots\E[x_n]^{m_n}
 \Big(1+\eta^{(m)}+\sum_{2\leq \sum\bar m_i< \sum m_i}\beta_{\bar m}\,\eta^{(\bar m)}\Big),
\end{align*}
where the summation is over normalized centered moments $\eta^{(\bar
  m)}$ of order two\footnote{The normalized centered moments of order
  one are always zero.} or larger and strictly smaller than $\sum m_i$;
and the $\beta_{\bar m}$ are appropriately selected nonnegative
constants.

\medskip

Conversely, we can also express the any normalized centered moments
$\eta^{(\bar m)}$ in terms of uncentered moment $\mu^{(\tilde m)}$ as
follows
\begin{align*}
  \eta^{(\bar m)}
  &=\E\Big[\Big(\frac{x_1-\E[x_1]}{\E[x_1]}\Big)^{\bar m_1}
  \Big(\frac{x_2-\E[x_2]}{\E[x_2]}\Big)^{\bar m_2} \cdots
  \Big(\frac{x_n-\E[x_n]}{\E[x_n]}\Big)^{\bar m_n}
  \Big]\\
  &=\E\Big[\Big(\frac{x_1}{\E[x_1]}-1\Big)^{\bar m_1}
  \Big(\frac{x_2}{\E[x_2]}-1\Big)^{\bar m_2} \cdots
  \Big(\frac{x_n}{\E[x_n]}-1\Big)^{\bar m_n}
  \Big]\\
  &=\frac{\mu^{(\bar m)}}{\E[x_1]^{\bar m_1}\E[x_2]^{\bar m_2}\cdots\E[x_n]^{\bar m_n}}
  + \sum_{\sum\tilde m_i< \sum m_i}\gamma_{\tilde m}\,\frac{\mu^{(\tilde m)}}{\E[x_1]^{\tilde m_1}\E[x_2]^{\tilde m_2}\cdots\E[x_n]^{\tilde m_n}},
\end{align*}
where the summation is over uncentered moments $\mu^{(\tilde m)}$ of
order strictly smaller than $\sum \bar m_i$ and the $\gamma_{\tilde m}$
are appropriately selected constants.

\medskip

The final formulas for low-dispersion moment closure can be obtained
directly from expansions of the non-normalized centered moments
\begin{align*}
  \E\Big[(x_1-\E[x_1])^{\bar m_1}(x_2-\E[x_2])^{\bar m_2} \cdots
  (x_n-\E[x_n])^{\bar m_n}\Big]
  &=\mu^{(\bar m)}+\sum_{\sum\tilde m_i< \sum m_i}\hat \gamma_{\tilde m}(\E[x_1],\E[x_2],\dots,\E[x_n])\mu^{(\tilde m)},
\end{align*}
because setting to zero some normalized centered $\eta^{(m)}$ is equivalent
to setting to zero the left-hand side of the equation above, from
which we conclude that
\begin{align*}
  \mu^{(\bar m)}\approx-\sum_{\sum\tilde m_i< \sum m_i} \hat \gamma_{\tilde m}(\E[x_1],\E[x_2],\dots,\E[x_n])\mu^{(\tilde m)}.
\end{align*}


\section{Quasi-deterministic Normalized Centered moments}\label{se:ap-quasi-deterministic}

The uncentered moment $\mu^{(m)}$ can be expressed in terms of the
quasi-deterministic normalized centered moment $\hat\eta^{(\bar m)}$ as follows
\begin{align*}
  \mu^{(m)} 
  &= \E[ x_1^{m_1} x_2^{m_2}\cdots x_n^{m_n} ] \\
  &= \phi_1^{m_1}\phi_2^{m_2}\cdots\phi_n^{m_n}
  \E\Big[\Big(1+\frac{x_1-\phi_1}{\phi_1}\Big)^{m_1}
  \Big(1+\frac{x_2-\phi_2}{\phi_2}\Big)^{m_2}\cdots
  \Big(1+\frac{x_n-\phi_n}{\phi_n}\Big)^{m_n}\Big]\\
  &=\phi_1^{m_1}\phi_2^{m_2}\cdots\phi_n^{m_n}
 \Big(1+\hat\eta^{(m)}+\sum_{2\leq \sum\bar m_i< \sum m_i}\beta_{\bar m}\,\hat\eta^{(\bar m)}\Big),
\end{align*}
where the summation is over quasi-deterministic normalized centered
moments $\hat\eta^{(\bar m)}$ of order one\footnote{Contrary to what happens
  for normalized centered moments, the quasi-deterministic normalized
  centered moments of order one are not necessarily zero.} or larger
and strictly smaller than $\sum m_i$; and the $\beta_{\bar m}$ are
appropriately selected nonnegative constants.

\medskip

Conversely, we can also express the any quasi-deterministic normalized
centered moments $\hat\eta^{(\bar m)}$ in terms of uncentered moment
$\mu^{(\tilde m)}$ as follows
\begin{align*}
  \hat\eta^{(\bar m)}
  &=\E\Big[\Big(\frac{x_1-\phi_1}{\phi_1}\Big)^{\bar m_1}
  \Big(\frac{x_2-\phi_2}{\phi_2}\Big)^{\bar m_2} \cdots
  \Big(\frac{x_n-\phi_n}{\phi_n}\Big)^{\bar m_n}
  \Big]\\
  &=\E\Big[\Big(\frac{x_1}{\phi_1}-1\Big)^{\bar m_1}
  \Big(\frac{x_2}{\phi_2}-1\Big)^{\bar m_2} \cdots
  \Big(\frac{x_n}{\phi_n}-1\Big)^{\bar m_n}
  \Big]\\
  &=\frac{\mu^{(\bar m)}}{\phi_1^{\bar m_1}\phi_2^{\bar m_2}\cdots\phi_n^{\bar m_n}}
  + \sum_{\sum\tilde m_i< \sum m_i}\gamma_{\tilde m}\,\frac{\mu^{(\tilde m)}}{\phi_1^{\tilde m_1}\phi_2^{\tilde m_2}\cdots\phi_n^{\tilde m_n}},
\end{align*}
where the summation is over uncentered moments $\mu^{(\tilde m)}$ of
order strictly smaller than $\sum \bar m_i$ and the $\gamma_{\tilde m}$
are appropriately selected constants.

\medskip

The final formulas for quasi-deterministic moment closure can be
obtained directly from expansions of the quasi-deterministic
non-normalized centered moments
\begin{align*}
  \E\Big[(x_1-\phi_1)^{\bar m_1}(x_2-\phi_2)^{\bar m_2} \cdots
  (x_n-\phi_n)^{\bar m_n}\Big]
  &=\mu^{(\bar m)}+\sum_{\sum\tilde m_i< \sum m_i}\hat \gamma_{\tilde m}(\phi_1,\phi_2,\dots,\phi_n)\mu^{(\tilde m)},
\end{align*}
because setting to zero some quasi-deterministic normalized centered
$\hat\eta^{(m)}$ is equivalent to setting to zero the left-hand side
of the equation above, from which we conclude that
\begin{align*}
  \mu^{(\bar m)}\approx-\sum_{\sum\tilde m_i< \sum m_i} \hat \gamma_{\tilde m}(\phi_1,\phi_2,\dots,\phi_n)\mu^{(\tilde m)}.
\end{align*}

%\bibliographystyle{ieeetr}
\bibliographystyle{abbrvnat}
\bibliography{strings,biology,software,mysoftware,jph,crossrefs}

\printindex

\end{document}
