update kinetic formalism and add competion start

[ool/lipid_simulation_formalism.git] / kinetic_formalism.Rnw

\documentclass[english,12pt]{article}
\usepackage{fancyhdr}
%\usepackage[pdftex]{graphicx}
\usepackage{graphicx}
\usepackage[bf]{caption2}
\usepackage{rotating}
\usepackage{multirow}
\usepackage{textcomp}
\usepackage{mathrsfs}
\usepackage{amssymb}
\usepackage{setspace}
\usepackage{txfonts}
\usepackage[light,all]{draftcopy}
\usepackage{fancyref}
\usepackage[hyperfigures,backref,bookmarks,colorlinks]{hyperref}
\usepackage[sectionbib,sort&compress,square,numbers]{natbib}
\usepackage[margin,inline,draft]{fixme}
\usepackage[x11names,svgnames]{xcolor}
\usepackage{texshade}
\newenvironment{narrow}[2]{%
  \begin{list}{}{%
      \setlength{\topsep}{0pt}%
      \setlength{\leftmargin}{#1}%
      \setlength{\rightmargin}{#2}%
      \setlength{\listparindent}{\parindent}%
      \setlength{\itemindent}{\parindent}%
      \setlength{\parsep}{\parskip}}%
  \item[]}{\end{list}}
\newenvironment{paperquote}{%
  \begin{quote}%
     \it
  }%
  {\end{quote}}
\renewcommand{\textfraction}{0.15}
\renewcommand{\topfraction}{0.85}
\renewcommand{\bottomfraction}{0.65}
\renewcommand{\floatpagefraction}{0.60}
%\renewcommand{\baselinestretch}{1.8}
\newenvironment{enumerate*}%
  {\begin{enumerate}%
    \setlength{\itemsep}{0pt}%
    \setlength{\parskip}{0pt}}%
  {\end{enumerate}}
\newenvironment{itemize*}%
  {\begin{itemize}%
    \setlength{\itemsep}{0pt}%
    \setlength{\parskip}{0pt}}%
  {\end{itemize}}
\oddsidemargin 0.0in 
\textwidth 6.5in
\raggedbottom
\clubpenalty = 10000
\widowpenalty = 10000
\pagestyle{fancy}
\author{Don Armstrong}
\title{OOL Kinetic Formalisms}
%\date{}
\onehalfspacing
\begin{document}
%\maketitle

<<results=hide,echo=FALSE>>=
require(lattice)
require(grid)
# R in cal / mol K
to.kcal <- function(k,temp=300) {
  gasconst <- 1.985
  return(-gasconst*temp*log(k)/1000)
}
@ 

\section{State Equation}
% double check this with the bits in the paper

Given a base forward kinetic parameter for the $i$th specie $k_{fi}$
(which is dependent on lipid type, that is PC, PE, PS, etc.), an
adjustment parameter $k_{fi\mathrm{adj}}$ based on the vesicle and the
specific specie (length, unsaturation, etc.) (see~\fref{eq:kf_adj}),
the molar concentration of monomer of the $i$th specie
$\left[C_{i_\mathrm{monomer}}\right]$, the surface area of the vesicle
$S_\mathrm{ves}$, the base backwards kinetic parameter for the $i$th
specie $k_{bi}$ which is also dependent on lipid type, its adjustment
parameter $k_{bi\mathrm{adj}}$ (see~\fref{eq:kb_adj}), and the molar
concentration of the $i$th specie in the vesicle
$\left[C_{i_\mathrm{ves}}\right]$, the change in concentration of the
$i$th specie in the vesicle per change in time $\frac{d
  C_{i_\mathrm{ves}}}{dt}$ can be calculated:

\begin{equation}
  \frac{d C_{i_\mathrm{ves}}}{dt} = k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]S_\mathrm{ves} -
  k_{bi}k_{bi\mathrm{adj}}C_{i_\mathrm{ves}}
  \label{eq:state}
\end{equation}

For $k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]$,
$k_{fi}$ has units of $\frac{\mathrm{m}}{\mathrm{s}}$,
$k_{fi\mathrm{adj}}$ and $k_{bi\mathrm{adj}}$ are unitless,
concentration is in units of $\frac{\mathrm{n}}{\mathrm{L}}$, surface
area is in units of $\mathrm{m}^2$, $k_{bi}$ has units of
$\frac{1}{\mathrm{s}}$ and $C_{i_\mathrm{ves}}$ has units of
$\mathrm{n}$, Thus, we have

\begin{equation}
  \frac{\mathrm{n}}{\mathrm{s}} = \frac{\mathrm{m}}{\mathrm{s}} \frac{\mathrm{n}}{\mathrm{L}} \mathrm{m}^2 \frac{1000\mathrm{L}}{\mathrm{m}^3} - 
  \frac{1}{\mathrm{s}} \mathrm{n}
  =
  \frac{\mathrm{m^3}}{\mathrm{s}} \frac{\mathrm{n}}{\mathrm{L}} \frac{1000\mathrm{L}}{\mathrm{m}^3} - \frac{\mathrm{n}}{\mathrm{s}}=
  \frac{\mathrm{n}}{\mathrm{s}} = 1000 \frac{\mathrm{n}}{\mathrm{s}} - \frac{\mathrm{n}}{\mathrm{s}}
  \label{eq:state_units}
\end{equation}

The 1000 isn't in \fref{eq:state} above, because it is unit-dependent.

\subsection{Forward adjustments ($k_{fi\mathrm{adj}}$)}
\label{sec:forward_adj}

The forward rate constant adjustment, $k_{fi\mathrm{adj}}$ takes into
account unsaturation ($un_f$), charge ($ch_f$), curvature ($cu_f$),
length ($l_f$), and complex formation ($CF1_f$), each of which are
modified depending on the specific specie and the vesicle into which
the specie is entering.

\begin{equation}
  k_{fi\mathrm{adj}} = un_f \cdot ch_f \cdot cu_f \cdot l_f \cdot CF1_f
  \label{eq:kf_adj}
\end{equation}

\newpage
\subsubsection{Unsaturation Forward}

In order for a lipid to be inserted into a membrane, a void has to be
formed for it to fill. Voids can be generated by the combination of
unsaturated and saturated lipids forming herterogeneous domains. Void
formation is increased when the unsaturation of lipids in the vesicle
is widely distributed; in other words, the insertion of lipids into
the membrane is greater when the standard deviation of the
unsaturation is larger. Assuming that an increase in width of the
distribution linearly decreases the free energy of activation, the
$un_f$ parameter must follow
$a^{\mathrm{stdev}\left(un_\mathrm{ves}\right)}$ where $a > 1$, so a
convenient starting base for $a$ is $2$:

\begin{equation}
  un_f = 2^{\mathrm{stdev}\left(un_\mathrm{ves}\right)}
  \label{eq:unsaturation_forward}
\end{equation}

The most common $\mathrm{stdev}\left(un_\mathrm{ves}\right)$ is around
$1.5$, which leads to a $\Delta \Delta G^\ddagger$ of
$\Sexpr{format(digits=3,to.kcal(2^1.5))}
\frac{\mathrm{kcal}}{\mathrm{mol}}$.

It is not clear that the unsaturation of the inserted monomer will
affect the rate of the insertion positively or negatively, so we do
not include a term for it in this formalism.


\setkeys{Gin}{width=3.2in}
<<fig=TRUE,echo=FALSE,results=hide,width=5,height=5>>=
curve(2^x,from=0,to=sd(c(0,4)),
      main="Unsaturation Forward",
      xlab="Standard Deviation of Unsaturation of Vesicle",
      ylab="Unsaturation Forward Adjustment")
@ 
<<fig=TRUE,echo=FALSE,results=hide,width=5,height=5>>=
curve(to.kcal(2^x),from=0,to=sd(c(0,4)),
      main="Unsaturation forward",
      xlab="Standard Deviation of Unsaturation of Vesicle",
      ylab="Unsaturation Forward (kcal/mol)")
@ 


\newpage
\subsubsection{Charge Forward}

A charged lipid such as PS approaching a vesicle with an average
charge of the same sign will experience repulsion, whereas those with
different signs will experience attraction, the degree of which is
dependent upon the charge of the monomer and the average charge of the
vesicle. If either the vesicle or the monomer has no charge, there
should be no effect of charge upon the rate. This leads us to the
following equation, $a^{-\left<ch_v\right> ch_m}$, where
$\left<ch_v\right>$ is the average charge of the vesicle, and $ch_m$
is the charge of the monomer. If either $\left<ch_v\right>$ or $ch_m$
is 0, the adjustment parameter is 1 (no change), whereas it decreases
if both are positive or negative, as the product of two real numbers
with the same sign is always positive. A convenient base for $a$ is
60, resulting in the following equation:


\begin{equation}
  ch_f = 60^{-\left<{ch}_v\right> {ch}_m}
  \label{eq:charge_forward}
\end{equation}

The most common $\left<{ch}_v\right>$ is around $-0.165$, which leads to
a range of $\Delta \Delta G^\ddagger$ from
$\Sexpr{format(digits=3,to.kcal(60^(-.165*-1)))}
\frac{\mathrm{kcal}}{\mathrm{mol}}$ to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$.

<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
x <- seq(-1,0,length.out=20)
y <- seq(-1,0,length.out=20)
grid <- expand.grid(x=x,y=y)
grid$z <- as.vector(60^(-outer(x,y)))
print(wireframe(z~x*y,grid,cuts=50,
                drape=TRUE,
                scales=list(arrows=FALSE),
                main="Charge Forward",
                xlab=list("Average Vesicle Charge",rot=30),
                ylab=list("Component Charge",rot=-35),
                zlab=list("Charge Forward",rot=93)))
rm(x,y,grid)
@ 
<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
x <- seq(-1,0,length.out=20)
y <- seq(-1,0,length.out=20)
grid <- expand.grid(x=x,y=y)
grid$z <- as.vector(to.kcal(60^(-outer(x,y))))
print(wireframe(z~x*y,grid,cuts=50,
                drape=TRUE,
                scales=list(arrows=FALSE),
                main="Charge Forward (kcal/mol)",
                xlab=list("Average Vesicle Charge",rot=30),
                ylab=list("Component Charge",rot=-35),
                zlab=list("Charge Forward (kcal/mol)",rot=93)))
rm(x,y,grid)
@ 


\newpage
\subsubsection{Curvature Forward}

Curvature is a measure of the intrinsic propensity of specific lipids
to form micelles (positive curvature), inverted micelles (negative
curvature), or planar sheets (zero curvature). In this formalism,
curvature is measured as the ratio of the size of the head to that of
the base, so negative curvature is bounded by $(0,1)$, zero curvature
is 1, and positive curvature is bounded by $(1,\infty)$. The curvature
can be transformed into the typical postive/negative mapping using
$\log$, which has the additional property of making the range of
positive and negative curvature equal, and distributed about 0.

As in the case of unsaturation, void formation is increased by the
presence of lipids with mismatched curvature. Thus, a larger
distribution of curvature in the vesicle increases the rate of lipid
insertion into the vesicle. However, a species with curvature $e^{-1}$
will cancel out a species with curvature $e$, so we have to log
transform (turning these into -1 and 1), then take the absolute value
(1 and 1), and finally measure the width of the distribution. Thus, by
using the log transform to make the range of the lipid curvature equal
between positive and negative, and taking the average to cancel out
exactly mismatched curvatures, we come to an equation with the shape
$a^{\left<\log cu_\mathrm{vesicle}\right>}$. A convenient base for $a$
is $10$, yielding:


\begin{equation}
 % cu_f = 10^{\mathrm{stdev}\left|\log cu_\mathrm{vesicle}\right|}
  cu_f = 10^{\left|\left<\log cu_\mathrm{vesicle} \right>\right|\mathrm{stdev} \left|\log cu_\mathrm{vesicle}\right|}
  \label{eq:curvature_forward}
\end{equation}

The most common $\left|\left<\log {cu}_v\right>\right|$ is around
$0.013$, which with the most common $\mathrm{stdev} \log
cu_\mathrm{vesicle}$ of $0.213$ leads to a $\Delta \Delta G^\ddagger$
of $\Sexpr{format(digits=3,to.kcal(10^(0.13*0.213)))}
\frac{\mathrm{kcal}}{\mathrm{mol}}$. This is a consequence of the
relatively matched curvatures in our environment.

% 1.5 to 0.75 3 to 0.33
<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
                      sd(abs(log(c(1,0.33)))),sd(abs(log(c(0.33,3)))))),length.out=20),
                    y=seq(0,max(c(mean(log(c(1,3)),
                      mean(log(c(1,0.33))),
                      mean(log(c(0.33,3)))))),length.out=20))
grid$z <- 10^(grid$x*grid$y)
print(wireframe(z~x*y,grid,cuts=50,
          drape=TRUE,
          scales=list(arrows=FALSE),
          xlab=list("Vesicle stdev log curvature",rot=30),
          ylab=list("Vesicle average log curvature",rot=-35),
          zlab=list("Vesicle Curvature Forward",rot=93)))
rm(grid)
@ 
<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
                      sd(abs(log(c(1,0.33)))),sd(abs(log(c(0.33,3)))))),length.out=20),
                    y=seq(0,max(c(mean(log(c(1,3)),
                      mean(log(c(1,0.33))),
                      mean(log(c(0.33,3)))))),length.out=20))
grid$z <- to.kcal(10^(grid$x*grid$y))
print(wireframe(z~x*y,grid,cuts=50,
          drape=TRUE,
          scales=list(arrows=FALSE),
          xlab=list("Vesicle stdev log curvature",rot=30),
          ylab=list("Vesicle average log curvature",rot=-35),
          zlab=list("Vesicle Curvature Forward (kcal/mol)",rot=93)))
rm(grid)
@ 

\newpage
\subsubsection{Length Forward}

As in the case of unsaturation, void formation is easier when vesicles
are made up of components of widely different lengths. Thus, when the
width of the distribution of lengths is larger, the forward rate
should be greater as well, leading us to an equation of the form
$x^{\mathrm{stdev} l_\mathrm{ves}}$, where $\mathrm{stdev}
l_\mathrm{ves}$ is the standard deviation of the length of the
components of the vesicle, which has a maximum possible value of 8 and
a minimum of 0 in this set of experiments. A convenient base for $x$
is 2, leading to:

\begin{equation}
  l_f = 2^{\mathrm{stdev} l_\mathrm{ves}}
  \label{eq:length_forward}
\end{equation}

The most common $\mathrm{stdev} l_\mathrm{ves}$ is around $3.4$, which leads to
a range of $\Delta \Delta G^\ddagger$ of
$\Sexpr{format(digits=3,to.kcal(2^(3.4)))}
\frac{\mathrm{kcal}}{\mathrm{mol}}$.

While it could be argued that increased length of the monomer could
affect the rate of insertion into the membrane, it's not clear whether
it would increase (by decreasing the number of available hydrogen
bonds, for example) or decrease (by increasing the time taken to fully
insert the acyl chain, for example) the rate of insertion or to what
degree, so we do not take it into account in this formalism.


<<fig=TRUE,echo=FALSE,results=hide,width=7,height=5>>=
curve(2^x,from=0,to=sd(c(12,24)),
      main="Length forward",
      xlab="Standard Deviation of Length of Vesicle",
      ylab="Length Forward Adjustment")
@ 
<<fig=TRUE,echo=FALSE,results=hide,width=7,height=5>>=
curve(to.kcal(2^x),from=0,to=sd(c(12,24)),
      main="Length forward",
      xlab="Standard Deviation of Length of Vesicle",
      ylab="Length Forward Adjustment (kcal/mol)")
@ 


\subsubsection{Complex Formation}
There is no contribution of complex formation to the forward reaction
rate in the current formalism.

\begin{equation}
  CF1_f=1
  \label{eq:complex_formation_forward}
\end{equation}

\subsection{Backward adjustments ($k_{bi\mathrm{adj}}$)}

Just as the forward rate constant adjustment $k_{fi\mathrm{adj}}$
does, the backwards rate constant adjustment $k_{bi\mathrm{adj}}$
takes into account unsaturation ($un_b$), charge ($ch_b$), curvature
($cu_b$), length ($l_b$), and complex formation ($CF1_b$), each of
which are modified depending on the specific specie and the vesicle
into which the specie is entering:


\begin{equation}
  k_{bi\mathrm{adj}} = un_b \cdot ch_b \cdot cu_b \cdot l_b \cdot CF1_b
  \label{eq:kb_adj}
\end{equation}

\subsubsection{Unsaturation Backward}

Unsaturation also influences the ability of a lipid molecule to leave
a membrane. If a molecule has an unsaturation level which is different
from the surrounding membrane, it will be more likely to leave the
membrane. The more different the unsaturation level is, the greater
the propensity for the lipid molecule to leave. However, a vesicle
with some unsaturation is more favorable for lipids with more
unsaturation than the equivalent amount of less unsatuturation, so the
difference in energy between unsaturation is not linear. Therefore, an
equation with the shape
$x^{\left| y^{-\left< un_\mathrm{ves}\right> }-y^{-un_\mathrm{monomer}} \right| }$
where $\left<un_\mathrm{ves}\right>$ is the average unsaturation of
the vesicle, and $un_\mathrm{monomer}$ is the average unsaturation. In
this equation, as the average unsaturation of the vesicle is larger,

\begin{equation}
  un_b = 7^{1-\left(20\left(2^{-\left<un_\mathrm{vesicle} \right>} - 2^{-un_\mathrm{monomer}} \right)^2+1\right)^{-1}}
  \label{eq:unsaturation_backward}
\end{equation}

The most common $\left<un_\mathrm{ves}\right>$ is around $1.7$, which leads to
a range of $\Delta \Delta G^\ddagger$ from
$\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-0)^2+1))))}
\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with 0 unsaturation
to
$\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-4)^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
for monomers with 4 unsaturations.


<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
grid <- expand.grid(x=seq(0,4,length.out=20),
                    y=seq(0,4,length.out=20))
grid$z <- (7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1)))
print(wireframe(z~x*y,grid,cuts=50,
          drape=TRUE,
          scales=list(arrows=FALSE),
          xlab=list("Average Vesicle Unsaturation",rot=30),
          ylab=list("Monomer Unsaturation",rot=-35),
          zlab=list("Unsaturation Backward",rot=93)))
rm(grid)
@ 
<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
grid <- expand.grid(x=seq(0,4,length.out=20),
                    y=seq(0,4,length.out=20))
grid$z <- to.kcal((7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1))))
print(wireframe(z~x*y,grid,cuts=50,
          drape=TRUE,
          scales=list(arrows=FALSE),
          xlab=list("Average Vesicle Unsaturation",rot=30),
          ylab=list("Monomer Unsaturation",rot=-35),
          zlab=list("Unsaturation Backward (kcal/mol)",rot=93)))
rm(grid)
@ 



\newpage
\subsubsection{Charge Backwards}
As in the case of monomers entering a vesicle, monomers leaving a
vesicle leave faster if their charge has the same sign as the average
charge vesicle. An equation of the form $ch_b = a^{\left<ch_v\right>
  ch_m}$ is then appropriate, and using a base of $a=20$ yields:

\begin{equation}
  ch_b = 20^{\left<{ch}_v\right> {ch}_m}
  \label{eq:charge_backwards}
\end{equation}

The most common $\left<ch_v\right>$ is around $-0.164$, which leads to
a range of $\Delta \Delta G^\ddagger$ from
$\Sexpr{format(digits=3,to.kcal(20^(-.164*-1)))}
\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with charge $-1$ to
$0\frac{\mathrm{kcal}}{\mathrm{mol}}$
for monomers with charge $0$.


<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
x <- seq(-1,0,length.out=20)
y <- seq(-1,0,length.out=20)
grid <- expand.grid(x=x,y=y)
grid$z <- as.vector(20^(outer(x,y)))
print(wireframe(z~x*y,grid,cuts=50,
          drape=TRUE,
          scales=list(arrows=FALSE),
          xlab=list("Average Vesicle Charge",rot=30),
          ylab=list("Component Charge",rot=-35),
          zlab=list("Charge Backwards",rot=93)))
rm(x,y,grid)
@ 
<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
x <- seq(-1,0,length.out=20)
y <- seq(-1,0,length.out=20)
grid <- expand.grid(x=x,y=y)
grid$z <- to.kcal(as.vector(20^(outer(x,y))))
print(wireframe(z~x*y,grid,cuts=50,
          drape=TRUE,
          scales=list(arrows=FALSE),
          xlab=list("Average Vesicle Charge",rot=30),
          ylab=list("Component Charge",rot=-35),
          zlab=list("Charge Backwards (kcal/mol)",rot=93)))
rm(x,y,grid)
@ 

\newpage
\subsubsection{Curvature Backwards}

The less a monomer's intrinsic curvature matches the average curvature
of the vesicle in which it is in, the greater its rate of efflux. If
the difference is 0, $cu_f$ needs to be one. To map negative and
positive curvature to the same range, we also need take the logarithm.
Increasing mismatches in curvature increase the rate of efflux, but
asymptotically.  An equation which satisfies this critera has the
form $cu_f = a^{1-\left(b\left( \left< \log cu_\mathrm{vesicle} \right>
      -\log cu_\mathrm{monomer}\right)^2+1\right)^{-1}}$. An
alternative form would use the aboslute value of the difference,
however, this yields a cusp and sharp increase about the curvature
equilibrium, which is decidedly non-elegant. We have chosen bases of
$a=7$ and $b=20$.

\begin{equation}
  cu_b = 7^{1-\left(20\left(\left< \log cu_\mathrm{vesicle} \right> -\log cu_\mathrm{monomer} \right)^2+1\right)^{-1}}
  \label{eq:curvature_backwards}
\end{equation}

The most common $\left<\log cu_\mathrm{ves}\right>$ is around $-0.013$, which leads to
a range of $\Delta \Delta G^\ddagger$ from
$\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(0.8))^2+1))))}
\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature 0.8
to
$\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(1.3))^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
for monomers with curvature 1.3 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature near 1.


<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
grid <- expand.grid(x=seq(0.8,1.33,length.out=20),
                    y=seq(0.8,1.33,length.out=20))
grid$z <- 7^(1-1/(20*(log(grid$x)-log(grid$y))^2+1))
print(wireframe(z~x*y,grid,cuts=50,
          drape=TRUE,
          scales=list(arrows=FALSE),
          xlab=list("Vesicle Curvature",rot=30),
          ylab=list("Monomer Curvature",rot=-35),
          zlab=list("Curvature Backward",rot=93)))
rm(grid)
@ 
<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
grid <- expand.grid(x=seq(0.8,1.33,length.out=20),
                    y=seq(0.8,1.33,length.out=20))
grid$z <- to.kcal(7^(1-1/(20*(log(grid$x)-log(grid$y))^2+1)))
print(wireframe(z~x*y,grid,cuts=50,
          drape=TRUE,
          scales=list(arrows=FALSE),
          xlab=list("Vesicle Curvature",rot=30),
          ylab=list("Monomer Curvature",rot=-35),
          zlab=list("Curvature Backward (kcal/mol)",rot=93)))
rm(grid)
@ 

\newpage
\subsubsection{Length Backwards}

In a model membrane, the dissociation constant increases by a factor
of approximately 3.2 per carbon decrease in acyl chain length (Nichols
1985). Unfortunatly, the experimental data known to us only measures
chain length less than or equal to the bulk lipid, and does not exceed
it, and is only known for one bulk lipid species (DOPC). We assume
then, that the increase is in relationship to the average vesicle, and
that lipids with larger acyl chain length will also show an increase
in their dissociation constant.

\begin{equation}
  l_b = 3.2^{\left|\left<l_\mathrm{ves}\right>-l_\mathrm{monomer}\right|}
  \label{eq:length_backward}
\end{equation}

The most common $\left<\log l_\mathrm{ves}\right>$ is around $17.75$, which leads to
a range of $\Delta \Delta G^\ddagger$ from
$\Sexpr{format(digits=3,to.kcal(3.2^abs(12-17.75)))}
\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with length 12
to
$\Sexpr{format(digits=3,to.kcal(3.2^abs(24-17.75)))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
for monomers with length 24 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature near 18.


<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
grid <- expand.grid(x=seq(12,24,length.out=20),
                    y=seq(12,24,length.out=20))
grid$z <- 3.2^(abs(grid$x-grid$y))
print(wireframe(z~x*y,grid,cuts=50,
          drape=TRUE,
          scales=list(arrows=FALSE),
          xlab=list("Average Vesicle Length",rot=30),
          ylab=list("Monomer Length",rot=-35),
          zlab=list("Length Backward",rot=93)))
rm(grid)
@ 
<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
grid <- expand.grid(x=seq(12,24,length.out=20),
                    y=seq(12,24,length.out=20))
grid$z <- to.kcal(3.2^(abs(grid$x-grid$y)))
print(wireframe(z~x*y,grid,cuts=50,
          drape=TRUE,
          scales=list(arrows=FALSE),
          xlab=list("Average Vesicle Length",rot=30),
          ylab=list("Monomer Length",rot=-35),
          zlab=list("Length Backward (kcal/mol)",rot=93)))
rm(grid)
@ 


\newpage
\subsubsection{Complex Formation Backward}

Complex formation describes the interaction between CHOL and PC or SM,
where PC or SM protects the hydroxyl group of CHOL from interactions
with water, the ``Umbrella Model''. PC ($CF1=2$) can interact with two
CHOL, and SM ($CF1=3$) with three CHOL ($CF1=-1$). If the average of
$CF1$ is positive (excess of SM and PC with regards to complex
formation), species with negative $CF1$ (CHOL) will be retained. If
average $CF1$ is negative, species with positive $CF1$ are retained.
An equation which has this property is
$CF1_b=a^{\left<CF1_\mathrm{ves}\right>
  CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{ves}\right>
    CF1_\mathrm{monomer}\right|}$, where difference of the exponent is
zero if the average $CF1$ and the $CF1$ of the specie have the same
sign, or double the product if the signs are different. A convenient
base for $a$ is $1.5$.


\begin{equation}
  CF1_b=1.5^{\left<CF1_\mathrm{ves}\right> CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{ves}\right> CF1_\mathrm{monomer}\right|}
  \label{eq:complex_formation_backward}
\end{equation}

The most common $\left<CF1_\mathrm{ves}\right>$ is around $0.925$,
which leads to a range of $\Delta \Delta G^\ddagger$ from
$\Sexpr{format(digits=3,to.kcal(1.5^(0.925*-1-abs(0.925*-1))))}
\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
formation $-1$ to
$\Sexpr{format(digits=3,to.kcal(1.5^(0.925*2-abs(0.925*2))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
for monomers with complex formation $2$ to
$0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
formation $0$.


<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
grid <- expand.grid(x=seq(-1,3,length.out=20),
                    y=seq(-1,3,length.out=20))
grid$z <- 1.5^(grid$x*grid$y-abs(grid$x*grid$y))
print(wireframe(z~x*y,grid,cuts=50,
          drape=TRUE,
          scales=list(arrows=FALSE),
          xlab=list("Vesicle Complex Formation",rot=30),
          ylab=list("Monomer Complex Formation",rot=-35),
          zlab=list("Complex Formation Backward",rot=93)))
rm(grid)
@ 
<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
grid <- expand.grid(x=seq(-1,3,length.out=20),
                    y=seq(-1,3,length.out=20))
grid$z <- to.kcal(1.5^(grid$x*grid$y-abs(grid$x*grid$y)))
print(wireframe(z~x*y,grid,cuts=50,
          drape=TRUE,
          scales=list(arrows=FALSE),
          xlab=list("Vesicle Complex Formation",rot=30),
          ylab=list("Monomer Complex Formation",rot=-35),
          zlab=list("Complex Formation Backward (kcal/mol)",rot=93)))
rm(grid)
@ 



\subsection{Per-Lipid Kinetic Parameters}

Each of the 5 lipid types have different kinetic parameters; to the
greatest extent possible, we have derived these from literature.

\begin{table}
  \centering
  \begin{tabular}{c c c c c c c}
    Type & $k_f$ & $k_b$ & Area (\r{A}$^2$) & Charge & CF1 & Curvature \\
    \hline
    PC   & $3.7\cdot 10^6$ & $2\cdot 10^{-5}$   & 63 & 0  & 2  & 0.8  \\
    PS   & $3.7\cdot 10^6$ & $1.5\cdot 10^{-5}$ & 54 & -1 & 0  & 1    \\
    CHOL & $5.1\cdot 10^7$ & $2.8\cdot 10^{-4}$ & 38 & 0  & -1 & 1.21 \\
    SM   & $3.7\cdot 10^6$ & $3.1\cdot 10^{-3}$ & 51 & 0  & 3  & 0.8  \\
    PE   & $2.3\cdot 10^6$ & $10^{-5}$          & 55 & 0  & 0  & 1.33 \\
  \end{tabular}
  \caption{Kinetic parameters of lipid types}
  \label{tab:kinetic_parameters_lipid_types}
\end{table}

\subsubsection{$k_f$ for lipid types}
For PC, $k_f$ was measured by Nichols85 to be $3.7\cdot 10^6
\frac{1}{\mathrm{M}\cdot \mathrm{s}}$ by the partitioning of
P-C$_6$-NBD-PC between DOPC vesicles and water. The method utilized by
Nichols85 has the weakness of using NBD-PC, with associated label
perturbations. As similar measures do not exist for SM or PS, we
assume that they have the same $k_f$. For CHOL, Estronca07 found a
value for $k_f$ of $5.1\cdot 10^7 \frac{1}{\mathrm{M}\cdot
  \mathrm{s}}$. For PE, Abreu04 found a value for $k_f$ of $2.3\cdot
10^6$. \fixme{I'm missing the notes on these last two papers, so this
  isn't correct yet.}

\subsubsection{$k_b$ for lipid types}

$k_b$ for PC was measured by Wimley90 using a radioactive label and
large unilammelar vesicles at 30\textdegree C. The other values were
calculated from the experiments of Nichols82 where the ratio of $k_b$
of different types was measured to that of PC.
See~\fref{tab:kinetic_parameters_lipid_types}.

assigned accordingly. kb(PS) was assumed to be the same as kb(PG) given
by Nichols82 (also ratio from kb(PC)).
kb(SM) is taken from kb(PC) of Wimley90 (radioactive), and then a ratio of
kb(PC)/kb(SM) taken from Bai97: = 34/2.2 = 15.45; 2.0 x 10-4 x 15.45 = 3.1 x
10-3 s -1.
kb(CHOL) taken from Jones90 (radioactive; POPC LUV; 37°).


\subsubsection{Area for lipid types}


\section{Simulation Methodology}

\subsection{Overall Architecture}

The simulation is currently run by single program written in perl
using various different modules to handle the subsidiary parts. It
produces output for each generation, including the step immediately
preceeding and immediately following a vesicle split (and optionally,
each step) that is written to a state file which contains the state of
the vesicle, environment, kinetic parameters, program invocation
options, time, and various other parameters necessary to recreate the
state vector at a given time. This output file is then read by a
separate program in perl to produce different output which is
generated out-of-band; output which includes graphs and statistical
analysis is performed using R (and various grid graphics modules)
called from the perl program.

The separation of simulation and output generation allows refining
output, and simulations performed on different versions of the
underlying code to be compared using the same analysis methods and
code. It also allows later simulations to be restarted from a specific
generation, utilizing the same environment.

Random variables of different distributions are calculated using
Math::Random, which is seeded for each run using entropy from the
linux kernel's urandom device.

\subsection{Environment Creation}

\subsubsection{Components}
The environment contains different concentrations of different
components. In the current set of simulations, there are
\Sexpr{1+4*5*7} different components, consisting of PC, PE, PS, SM,
and CHOL, with all lipids except for CHOL having 5 possible
unsaturations rangiong from 0 to 4, and 7 lengths from $12,14,...,22$
($7\cdot 5\cdot4+1=\Sexpr{1+4*5*7}$). In cases where the environment
has less than the maximum number of components, the components are
selected in order without replacement from a randomly shuffled deck of
component (with the exception of CHOL) represented once until the
desired number of unique components are obtained. CHOL is over
representated $\Sexpr{5*7}$ times to be at the level of other lipid
types, ensures that the probability of CHOL being asbent in the
environment is the same as the probability of one of the other lipid
types (PS, PE, etc.) being entirely absent. This reduces the number of
simulations with a small number of components which are entirely
devoid of CHOL.

\subsubsection{Concentration}
Once the components of the environment have been selected, the
concentration of thoes components is determined. In experiments where
the environmental concentration is the same across all lipid
components, the concentration is set to $10^{-10}\mathrm{M}$. In other
cases, the environmental concentration is set to a random number from
a gamma distribution with shape parameter 2 and an average of
$10^{-10}\mathrm{M}$. The base concentration ($10^{-10}\mathrm{M}$)
can also be altered in the initialization of the experiment to
specific values for specific lipid types or components.

\subsection{Initial Vesicle Creation}

Initial vesicles are seeded by selecting lipid molecules from the
environment until the vesicle reaches a specific starting size. The
vesicle starting size has gamma distribution with shape parameter 2
and a mean of the experimentally specified starting size, with a
minimum of 5 lipid molecules. Lipid molecules are then selected to be
added to the lipid membrane according to four different methods. In
the constant method, molecules are added in direct proportion to their
concentration in the environment. The uniform method adds molecules in
proportion to their concentration in the environment scaled by a
uniform random value, whereas the random method adds molecules in
proportion to a uniform random value. The final method is a binomial
method, which adjust the porbability of adding a molecule of a
specific component by the concentration of that component, and then
adds the molecules one by one to the membrane. This final method is
also used in the first three methods to add any missing molecules to
the starting vesicle which were unallocated due to the requirement to
add integer numbers of molecules.

\subsection{Simulation Step}

Once the environment has been created and the initial vesicle has been
formed, molecules join and leave the vesicle based on the kinetic
parameters and state equation discussed until the vesicle splits
forming two daughter vesicles, one of which the program continues to
follow.

\subsubsection{Calculation of Vesicle Properties}

\label{sec:ves_prop}
$S_\mathrm{ves}$ is the surface area of the vesicle, and is the sum of
the surface area of all of the individual lipid components; each lipid
type has a different surface area; we na\"ively assume that the lipid
packing is optimal, and there is no empty space.

\subsubsection{Joining and Leaving of Lipid Molecules}

Determining the number of molecules to add to the lipid membrane
($n_i$) requires knowing $k_{fi\mathrm{adj}}$, the surface area of the
vesicle $S_\mathrm{ves}$ (see~\fref{sec:ves_prop}), the time interval
$dt$ during which lipids are added, the base $k_{fi}$, and the
concentration of the monomer in the environment
$\left[C_{i\mathrm{ves}}\right]$ (see~\fref{eq:state}).
$k_{fi\mathrm{adj}}$ is calculated (see~\fref{eq:kf_adj}) based on the
vesicle properties and their hypothesized effect on the rate (in as
many cases as possible, experimentally based)
(see~\fref{sec:forward_adj} for details). $dt$ can be varied
(see~\fref{sec:step_duration}), but for a given step is constant. This
leads to the following:

$n_i = k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]S_\mathrm{ves}\mathrm{N_A}dt$

In the cases where $n_i > 1$, the integer number of molecules is
added. Fractional $n_i$ or the fractional remainder after the addition
of the integer molecules are the probability of adding a specific
molecule, and are compared to a uniformly distributed random value
between 0 and 1. If the random value is less than or equal to the
fractional part of $n_i$, an additional molecule is added.

Molecules leaving the vesicle are handled in a similar manner, with 

$n_i = k_{bi}k_{bi\mathrm{adj}}C_{i_\mathrm{ves}}\mathrm{N_A}dt$.

While programatically, the molecule removal happens after the
addition, the properties that each operates on are the same, so they
can be considered to have been added and removed at the same instant.
[This also avoids cases where a removal would have resulted in
negative molecules for a particular type, which is obviously
disallowed.]

\subsubsection{Step duration}
\label{sec:step_duration}

$dt$ is the time taken for each step of joining, leaving, and checking
split. In the current implementation, it starts with a value of
$10^{-6}\mathrm{s}$ but this is modified in between each step if the
number of molecules joining or leaving is too large or small. If more
than half of the starting vesicle size molecules join or leave in a
single step, $dt$ is reduced by half. If less than the starting
vesicle size molecules join or leave in 100 steps, $dt$ is doubled.
This is necessary to curtail run times and to automatically adjust to
different experimental runs.

(In every run seen so far, the initial $dt$ was too small, and was
increased before the first generation occured; at no time was $dt$ too
large.)

\subsubsection{Vesicle splitting}

If a vesicle has grown to a size which is more than double the
starting vesicle size, the vesicle splits. More elaborate mechanisms
for determining whether a vesicle should split are of course possible,
but not currently implemented. Molecules are assigned to the daughter
vesicles at random, with each daughter vesicle having an equal chance
of getting a single molecule. The number of molecules to assign to the
first vesicle has binomial distribution with a probability of an event
in each trial of 0.5 and a number of trials equal to the number of
molecules.

\subsection{Output}

The environment, initial vesicle, and the state of the vesicle
immediately before and immediately after splitting are stored to disk
to produce later output.

\section{Analyzing output}

Analyzing of output is handled by a separate perl program which shares
many common modules with the simulation program. Current output
includes simulation progress, summary tables, summary statistics, and
various graphs.

\subsection{PCA plots}

Vesicles have many different axes which contribute to their variation
between subsequent generations; two major groups of axes are the
components and properties of vesicles. Each component in a vesicle is
an axis on its own; it can be measured either as an absolute number of
molecules in each component, or the fraction of molecules of that
component over the total number of molecules; the second approach is
often more convenient, as it allows vesicles of different number of
molecules to be more directly compared (though it hides any affect of
vesicle size).

In order to visualize the transition of subsequent generations of
vesicles from their initial state in the simulation, to their final
state at the termination of 

\subsection{Carpet plots}



% \bibliographystyle{plainnat}
% \bibliography{references.bib}


\end{document}