use knitr caching

[ool/lipid_simulation_formalism.git] / kinetic_formalism_competition.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{acronym}
\usepackage{array}
\usepackage{dcolumn}
\usepackage{booktabs}
\usepackage[noblocks]{authblk}
\usepackage[hyperfigures,bookmarks,colorlinks,citecolor=black,filecolor=black,linkcolor=black,urlcolor=black]{hyperref}
\usepackage[sectionbib,sort&compress,numbers]{natbib}
\usepackage[nomargin,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}}

\newcommand{\DLA}[1]{\textcolor{red}{\fxnote{DLA: #1}}}

\newcommand{\RZ}[1]{\textcolor{red}{\fxnote{RZ: #1}}}

\newcommand{\OM}[1]{\textcolor{red}{\fxnote{OM: #1}}}


%\SweaveOpts{prefix.string=pub_318_phys_bio_submission_figures_suppl/pub_318_phys_bio_sub_suppl_fig,pdf=TRUE,eps=TRUE}
\oddsidemargin 0.0in 
\textwidth 6.5in
\raggedbottom
\clubpenalty = 10000
\widowpenalty = 10000
\pagestyle{fancy}
\author{}
\title{Supplemental Material}
\date{}
% rubber: module bibtex
% rubber: module natbib
\onehalfspacing
\begin{document}
\maketitle

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

% \setcounter{figure}{0} \setcounter{table}{0} 
\makeatletter
\renewcommand{\thefigure}{S\@arabic\c@figure}
\renewcommand{\thetable}{S\@arabic\c@table}
\makeatother

\section{Competition Implementation}
\subsection{Implementation changes}

\begin{itemize}
\item settable maximum number of vesicles to track (default $10^4$)
\item start with 1~L ($10^{-3}$~m$^3$) cube
\item if at any point the number of vesicles exceeds the maximum
  number, chop the volume and environment molecule number into tenths,
  randomly select one tenth of the vesicles, and continue tracking.
\item generations will be counted per vesicle, and each progeny
  vesicle will have a generation number one greater than the parental
  vesicle.
\item 100 generations can result in as many as $2^{100}$
  ($\Sexpr{to.latex(format(digits=3,2^100))}$) vesicles or as few as
  101 vesicles.
\item Environment will use a specific number of each component instead
  of a constant concentration; as the number may be larger than
  \texttt{long long} ($2^{64}$), we use libgmp to handle an arbitrary
  precision number of components
\end{itemize}

\subsection{Infrastructure changes}
\begin{itemize}
\item Rewrite core bits in C
\item Use libgmp for handling large ints
\item Use openmpi to split the calculations out over multiple
  machines/processors and allow deploying to larger
  clusters/supercomputers
\end{itemize}



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

For a system with monomers $(_\mathrm{monomer})$ and a vesicle
$(_\mathrm{vesicle})$, the change in composition of the $i$th component of
a lipid vesicle per change in time ($d \left[C_{i_\mathrm{vesicle}}\right]/dt$)
can be described by a modification of the basic mass action law:

\begin{equation}
  \frac{d \left[C_{i_\mathrm{vesicle}}\right]}{dt} = k_{\mathrm{f}i}k_{\mathrm{f}i\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]\left[S_\mathrm{vesicle}\right] -
  k_{\mathrm{b}i}k_{\mathrm{b}i\mathrm{adj}}\left[C_{i_\mathrm{vesicle}}\right]
  \label{eq:state}
\end{equation}

The base forward kinetic parameter for the $i$th component is $k_{\mathrm{f}i}$
and is dependent on the particular lipid type (PC, PS, SM, etc.). The
forward adjustment parameter, $k_{\mathrm{f}i\mathrm{adj}}$, is based on the
properties of the vesicle and the specific component (type, length,
unsaturation, etc.) (see Equation~\ref{eq:kf_adj}, and
Section~\ref{sec:kinetic_adjustments}).
$\left[C_{i_\mathrm{monomer}}\right]$ is the molar concentration of
monomer of the $i$th component. $\left[S_\mathrm{vesicle}\right]$ is the surface
area of the vesicle per volume. The base backwards kinetic parameter
for the $i$th component is $k_{\mathrm{b}i}$ and its adjustment parameter
$k_{\mathrm{b}i\mathrm{adj}}$ (see Equation~\ref{eq:kb_adj}, and
Section~\ref{sec:kinetic_adjustments}).
$\left[C_{i_\mathrm{vesicle}}\right]$ is the molar concentration of
the $i$th component in the vesicle.

\subsection{Per-Lipid Kinetic Parameters}

<<echo=FALSE,results='hide'>>=
kf.prime <- c(3.7e6,3.7e6,5.1e7,3.7e6,2.3e6)
kf <- (as.numeric(kf.prime)*10^-3)/(63e-20*6.022e23)
@ 

Each of the 5 lipid types has different kinetic parameters; where
available, these were taken from literature.

% \begin{table}
%   \centering
%   \begin{tabular}{c c c c c c c c}
%     \toprule
%     Type & $k_\mathrm{f}$ $\left(\frac{\mathrm{m}}{\mathrm{s}}\right)$ & $k'_\mathrm{f}$ $\left(\frac{1}{\mathrm{M} \mathrm{s}}\right)$ & $k_\mathrm{b}$ $\left(\mathrm{s}^{-1}\right)$ & Area $\left({\AA}^2\right)$ & Charge & $\mathrm{CF}1$ & Curvature \\
%         \midrule
%     PC   & $\Sexpr{to.latex(format(digits=3,scientific=TRUE,kf[1]))}$ & $3.7 \times 10^6$ & $2   \times 10^{-5}$ & 63 & 0  & 2  & 0.8  \\
%     PS   & $\Sexpr{to.latex(format(digits=3,scientific=TRUE,kf[2]))}$ & $3.7 \times 10^6$ & $1.25\times 10^{-5}$ & 54 & -1 & 0  & 1    \\
%     CHOL & $\Sexpr{to.latex(format(digits=3,scientific=TRUE,kf[3]))}$ & $5.1 \times 10^7$ & $2.8 \times 10^{-4}$ & 38 & 0  & -1 & 1.21 \\
%     SM   & $\Sexpr{to.latex(format(digits=3,scientific=TRUE,kf[4]))}$ & $3.7 \times 10^6$ & $3.1 \times 10^{-3}$ & 61 & 0  & 3  & 0.8  \\
%     PE   & $\Sexpr{to.latex(format(digits=3,scientific=TRUE,kf[5]))}$ & $2.3 \times 10^6$ & $1   \times 10^{-5}$ & 55 & 0  & 0  & 1.33 \\
%     \bottomrule
%   \end{tabular}
%   \caption{Kinetic parameters and molecular properties of lipid types}
%   \label{tab:kinetic_parameters_lipid_types}
% \end{table}

%%% \DLA{I think we may just reduce these three sections; area, $k_\mathrm{f}$
%%%   and $k_\mathrm{b}$ to Table~\ref{tab:kinetic_parameters_lipid_types} with
%%%   references.}
%%% 
%%% \RZ{Yes, but we also have to have then as comments the numbers that
%%%   are not supported by refs}

\subsubsection{$k_\mathrm{f}$ for lipid types}
$k_{\mathrm{f}_\mathrm{PC}}$ was measured by
\citet{Nichols1985:phospholipid_monomer_vesicle_thermodynamics} and
was found to be $3.7\times 10^6$~$\frac{1}{\mathrm{M}\mathrm{s}}$ by the
partitioning of P-C$_6$-NBD-PC between DOPC vesicles and water. As
similar references are not available for SM or PS, we assume that they have
the same $k_\mathrm{f}$. For CHOL, no direct measurement of $k_\mathrm{f}$ is available,
however, \citet{Estronca2007:dhe_kinetics} measured the transfer of
DHE from BSA to LUV, and found a $k_\mathrm{f}$ of $5.1\times 10^7$~%
$\frac{1}{\mathrm{M} \mathrm{s}}$. We assume that this value is close
to that of CHOL, and use it for $k_{\mathrm{f}_\mathrm{CHOL}}$. In the case of
PE, \citet{Abreu2004:kinetics_ld_lo} measured the association of
NBD-DMPE with POPC LUV found a value for $k_\mathrm{f}$ of $2.3 \times 10^{6}$~%
$\frac{1}{\mathrm{M} \mathrm{s}}$. These three authors used a slightly
different kinetic formalism ($\frac{d\left[A_v\right]}{dt} =
k'_\mathrm{f}[A_m][L_v] - k_\mathrm{b}[A_v]$), so we correct their values of $k_\mathrm{f}$ by
multiplying by the surface area of a mole of lipids.

\subsubsection{$k_\mathrm{b}$ for lipid types}

\citet{Wimley1990:dmpc_exchange} measured the half time for the
exchange of [$^3$H]DMPC between LUV at 30\textdegree C, and found it
to be 9.6 hr. As this is a first order reaction, and the primary
limiting step in exchange is lipid desorption, $k_\mathrm{b}$ for DMPC is
$k_{\mathrm{b}_\mathrm{PC}}=\frac{\log 2}{9.6 \times 60 \times 60} \approx
\Sexpr{to.latex(format(log(2)/(9.6*60*60),digits=2))}
\mathrm{s}^{-1}$. We assume that $k_\mathrm{b}$ for SM is the same as for PC.
To estimate the $k_\mathrm{b}$ of PE and PS, we used the data from
\citet{Nichols1982:ret_amphiphile_transfer} who measured the rate of
exchange of the fluorescent label C$_6$-NBD attached to different
lipid species. Although the values of kb are different for the labeled
and unlabeled lipids, we assume that the ratios of the kinetics
constants for the lipid types are the same. Furthermore we assume that
PG behaves similarly to PS. Thus, we can determine the $k_\mathrm{b}$ of PE and
PS from the already known $k_\mathrm{b}$ of PC. For C$_6$-NBD labeled PC,
\citet{Nichols1982:ret_amphiphile_transfer} obtained a $k_\mathrm{b}$ of
$0.89$~$\mathrm{min}^{-1}$, PE of $0.45$~$\mathrm{min}^{-1}$ and PG of
$0.55$~$\mathrm{min}^{-1}$. Assuming the ratios of the rate of
exchange are the same for unlabeled lipids and labeled lipids, we can
determine the $k_\mathrm{b}$ of PE and PS from the already known $k_\mathrm{b}$ of
PC~\citep{Wimley1990:dmpc_exchange}. Calculating the ratio, this leads
us to $k_{\mathrm{b}_\mathrm{PE}} =
\frac{k_{\mathrm{b}_\mathrm{PC}}\times\mathrm{PE}}{\mathrm{PC}} \approx
\frac{2\times 10^{-5}\,\mathrm{s}^{-1} \times
  0.45\,\mathrm{min}^{-1}}{0.89\,\mathrm{min}^{-1}} \approx
\Sexpr{to.latex(format(log(2)/(9.6*60*60)*0.45/0.89,digits=2))}$~$\mathrm{s}^{-1}$
and likewise, $k_{\mathrm{b}_\mathrm{PS}}\approx
\Sexpr{to.latex(format(log(2)/(9.6*60*60)*0.55/0.89,digits=3))}$~$\mathrm{s}^{-1}$.
The $k_\mathrm{b}$ of SM was determined using the work of
\citet{Bai1997:lipid_movementbodipy}, who measured spontaneous
transfer of C$_5$-DMB-SM and C$_5$-DMB-PC from donor and acceptor
vesicles, finding $3.4\times10^{-2}$~$\mathrm{s}^{-1}$ and
$2.2\times10^{-3}$~$\mathrm{s}^{-1}$ respectively; using the ratio of
$k_\mathrm{b}$ of C$_5$-DMB-SM to the $k_\mathrm{b}$ of C$_5$-DMB-PC times the $k_\mathrm{b}$ of
PC ($\frac{3.4 \times 10^{-2}\mathrm{s}^{-1}}{2.2 \times
  10^{-3}\mathrm{s}^{-1}}
\Sexpr{to.latex(format(log(2)/(9.6*60*60),digits=2))}\mathrm{s}^{-1}$,
we obtain $k_{\mathrm{b}_\mathrm{SM}} \approx
\Sexpr{to.latex(format(log(2)/(9.6*60*60)*3.4e-2/2.2e-3,digits=2,scientific=TRUE))}$.
In the case of CHOL, \citet{Jones1990:lipid_transfer} measured the
$t_{1/2}$ of [$^3$H] transfer from POPC vesicles and found it to be 41
minutes, leading to a $k_{\mathrm{b}_\mathrm{CHOL}} = \frac{\log 2}{41\times
  60} \approx
\Sexpr{to.latex(format(log(2)/(41*60),digits=2,scientific=TRUE))}$~$\mathrm{s}^{-1}$.

\subsubsection{Headgroup Surface Area for lipid types}

% Smondyrev, A. M., and M. L. Berkowitz. 1999b. United atom force
% field for phospholipid membranes: constant pressure molecular
% dynamics simulation of dipalmitoylphosphatidicholine/water system.
% J. Comput. Chem. 50:531–545.

Different lipids have different headgroup surface areas, which contributes to
$\left[S_\mathrm{vesicle}\right]$. \citet{Smaby1997:pc_area_with_chol}
measured the surface area of POPC with a Langmuir film balance, and
found it to be 63~$\AA^2$ at $30$~$\frac{\mathrm{mN}}{\mathrm{m}}$.
Molecular dynamic simulations found an area of 54 $\AA^2$ for
DPPS\citep{Cascales1996:mds_dpps_area,Pandit2002:mds_dpps}, which is
in agreement with the experimental value of 56~$\AA^2$ found using a
Langmuir balance by \citet{Demel1987:ps_area}.
\citet{Shaikh2002:pe_phase_sm_area} measured the area of SM using a
Langmuir film balance, and found it to be 61~$\AA^2$. Using $^2$H NMR,
\citet{Thurmond1991:area_of_pc_pe_2hnmr} found the area of
DPPE-d$_{62}$ to be 55.4 $\AA^2$. \citet{Robinson1995:mds_chol_area}
found an area for CHOL of 38~$\AA^2$ using molecular dynamic
simulations.

% robinson's chol area is kind of crappy; they did it using MDS, but
% vaguely; they don't indicate what it was at the end, and they only
% had 2 molecules of CHOL.



% PC: 63 A (Smaby97rd) (Pandit (?)
% PS: 54 A (Pandit02LIB) (Cascales 1996; J. Chem. Phys 104:2713, Mukhopadhyay2004)
% CHOL:  38 A (Robinson 1995) (Lundberg 1982)
% SM: 51 A (Shaikh2002, but 61(?))
% PE: 55 A (Thurmond, 1991) (Pandit2002)
% 
% Compare to results by Petrache2004 with DOPC of 72.5, DOPS of 65.3.
% 
% The cross‑sectional area of cholesterol is ~37 ֵ2 (Lundberg, B. 1982. A
% surface film study of the lateral packing of phosphatidylcholine and
% cholesterol. Chem. Phys. Lipids. 31:23-32).
% 
% From Pandit02LIB: the area per headgroup for PE is ~10‑20\% (Thurmond
% et al., 1991) smaller than the area per PC molecule. For DLPE, it is
% 50.6 ֵ2 (McIntosh and Simon, 1986; Damodaran et al., 1992). The
% estimated area per molecule for DPPE is ~55.4 ֵ2 (Thurmond et al.,
% 1991).
% 
% From Pandit02LIB: (from their simulations) average area per headgroup
% for DPPS molecules is 53.75 ± 0.10 ֵ2. The values inferred from the
% experiments are in the region between 45 and 55 ֵ2 (Cevc et al., 1981.
% Titration of the phase transition of phosphatidylserine bilayer
% membranes. Effect of pH, surface electrostatics, ion binding, and
% head‑group hydration. Biochemistry. 20:4955‑4965; | Demel et al.,
% 1987. Monolayer characteristics and thermal behavior of natural and
% synthetic phosphatidylserines. Biochemistry. 26:8659‑8665). Note that
% the average area per headgroup for DPPC bilayer obtained from
% simulations and measurements is ~62 ֵ A2. One would expect that the
% area per headgroup in case of DPPS molecules will be larger than the
% area per DPPC molecule, because the headgroups of DPPS are charged and
% therefore repel each other. Contrary to this expectation, the area per
% headgroup in DPPS is ~13\% smaller than that of DPPC. López Cascales
% et al. (1996. Molecular dynamics simulation of a charged biological
% membrane. J. Chem. Phys. 104:2713‑2720.) proposed that this reduction
% in the area per headgroup is due to a strong intermolecular
% coordination between DPPS molecules. MD gave area per POPS of 55 A2 at
% 300K (Mukhopadhyay04LIB).
% 
% However, by 2H NMR and X-ray, Petrache04LIB determined the area of
% DOPS at 30 C to be 65.3 A2, considerably smaller than DOPC, determined
% by the same group in another paper to be 72.5 ֵ A2
% (Tristram‑Nagle00LIB).
% 
% From Shaikh02: At 35 mN/m of surface pressure, the areas/molecule for
% POPE, SM, and CHOL were found to be 65.2, 61.1, and 38.8 Å2.
% 
% In Lonchin99: POPC area 72 A2 ((a) Cornell, B. A.; Middlehurst, J.;
% Separovic, F. Biochim. Biophys. Acta1980, 598, 405−410), oleic acid
% area 32 A2 ((a) Small, D. M. In Handbook of Lipid Research; Hanahan,
% D. J., Ed.; Plenum Press:  New York, 1986; Vol. 4, Chapter 9.
% 
% (b) Jain, M. K. Introduction to Biological Membranes, 2nd ed.; John
% Wiley \& Sons:  New York, 1988; Chapter 3). Cistola, D. P.; Atkinson,
% D.; Hamilton, J. A.; Small, D. M. Biochemistry1986, 25, 2804−2812.
% 
% From Marsh96: The values obtained for the headgroup area per molecule,
% fitted as free parameters in Fig. 6, are APE = 58 A2 and APC = 78 A2,
% for DOPE and DOPC, respectively. These values are within the range
% obtained for the corresponding areas per molecule of
% phosphatidylethanolamines and phosphatidylcholines in the fluid
% lamellar phase (see, e.g., Marsh, 1990), and correspond to values for
% the lipid packing parameter of V/Al = 1.08 and 1.34 for DOPC and DOPE,
% respectively. The former value is consistent with DOPC alone forming
% hydrated lamellar structures rather than nonlamellar structures. The
% latter value lies at the lower end of the packing parameters obtained
% for phosphatidylethanolamines.
% 
% From Kumar91: The head group area chosen was 71.7 A2 for PC and
% lysoPC, 42 A2 for PE, and 19 A2 for Chol (5, 14). The value selected
% for PC is in agreement with the 71 A2 value determined for PCs with
% saturated acyl chains (14). The hydrocarbon volume and the length
% chosen for Chol are 400 A3 and 17.25 A, respectively (5).
% 
% 
% 
% From Sampaio05: Besides this work and our own earlier report on the
% association of NBD-DMPE with lipid bilayers (Abreu et al., 2004), we
% are aware of only one other report in the literature (Nichols, 1985)
% in which all the kinetic constants of lipid-derived amphiphile
% association with lipid bilayer membranes were experimentally measured.
% {indeed; everything is k- !!!; rz}
% 
% From McLean84LIB: Although it is difficult to measure cmc values for
% the sparingly soluble lipids used in this study, estimates for
% lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds,
% 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution;
% Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6
% X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X
% 10-8 M for DMPC was estimated from the linear relationship between ln
% cmc and the number of carbons in the PC acyl chain by using data for n
% = 7, 8, 10, and 16 [summarized in Tanford (1980)].
% 
% From Toyota08: Recently, several research groups have reported
% self-reproducing systems of giant vesicles that undergo a series of
% sequential transformations: autonomous growth, self-division, and
% chemical reactions to produce membrane constituents within the giant
% vesicles.44-47
% 
% Vesicle sizes of 25 nm for SUV and 150 nm for LUV were mentioned by
% Thomas02.
% 
% From Lund-Katz88: Charged and neutral small unilamellar vesicles
% composed of either saturated PC, unsaturated PC, or SM had similar
% size distributions with diameters of 23 \& 2 nm.
% 
% From Sampaio05LIB: The exchange of lipids and lipid derivatives
% between lipid bilayer vesicles has been studied for at least the last
% 30 years. Most of this work has examined the exchange of amphiphilic
% molecules between a donor and an acceptor population. The measured
% efflux rates were shown in almost all cases, not surprisingly, to be
% first order processes. In all of this work, the insertion rate has
% been assumed to be much faster than the efflux rate. Having measured
% both the insertion and desorption rate constants for amphiphile
% association with membranes, our results show that this assumption is
% valid. In several cases reported in the literature, the insertion rate
% constant was assumed, although never demonstrated, to be a
% diffusion-controlled process.
% 
% (for methods? From McLean84LIB: The activation free energies and free
% energies of transfer from self-micelles to water increase by 2.2 and
% 2.1 kJ mol-' per methylene group, respectively. {see if we can use it
%   to justify arranging our changed in activating energy around 1
%   kcal/mol; rz}).
% 
% Jones90 give diameter of LUV as 100 nm, and of SUV as 20 nm; that
% would give the number of molecules per outer leaflet of a vesicle as
% 1500.
% 
% Form Simard08: Parallel studies with SUV and LUV revealed that
% although membrane curvature does have a small effect on the absolute
% rates of FA transfer between vesicles, the ΔG of membrane desorption
% is unchanged, suggesting that the physical chemical properties which
% govern FA desorption are dependent on the dissociating molecule rather
% than on membrane curvature. However, disagreements on this fundamental
% issue continue (57, 61, 63, 64)
% 
% (methods regarding the curvature effect: Kleinfeld93 showed that the
% transfer parameters of long-chain FFA between the lipid vesicles
% depend on vesicle curvature and composition. Transfer of stearic acid
% is much slower from LUV as compared to SUV).
% 
% From McLean84: In a well-defined experimental system consisting of
% unilamellar lipid vesicles, in the absence of protein, the
% rate-limiting step for the overall exchange process is desorption
% (McLean \& Phillips, 1981). {thus I can take exchange data for the
%   estimation of k- rz; 8/11/08}.

\subsubsection{Complex Formation 1}

\citet{Thomas1988:chol_transfer} found that SM significantly decreases
the rate of cholesterol transfer from PC, while PE and PS have no
effect at physiologically significant levels. We attribute this to the
formation of complexes between SM and CHOL, which retards the rate of
efflux of CHOL from the membrane. Similar complexes exist between PC
and CHOL, where it was shown that two CHOL molecules complex with one
PC~\citep{Huang1999:chol_solubility_model,
  Huang1999:cholesterol_solubility,McConnell2006,McConnell2003}. It
has also been shown that SM binds more CHOL molecules than
PC~\citep{Katz1988:pl_packing_chol}. We assume that three CHOL complex
with one SM, leading to $\mathrm{CF}1$ values of $3$, $2$, and $-1$
for SM, PC, and CHOL, respectively.

\subsubsection{Curvature}

We used the formulation for curvature of
\citet{Israelachvili1975:amphiphile_self_assembly}, or $S=\frac{v}{l_c
  a}$, where $S$ is the curvature, which ranges from $0$ to $\infty$,
$l_c$ is the critical length of the acyl chain, $v$ is the volume of
the lipid, and $a$ is the area of the head group.
\citet{Kumar1991:lipid_packing} found a value for $S$ of 1.33 for PE,
1.21 for CHOL, and 0.8 for diC$_{16}$ PC. We assume that PS has neutral
curvature (value of 1), and that SM has the same curvature as PC
(0.8). It should also be noted that in reality the curvature parameter
changes with length, but at longer chain lengths, is effectively
constant; in the current model, curvature is held constant for each
species.

% 
% 1. Israelachvili, J. N., Mitchell, D. J. & Ninham, B. W. (1976) J.
%    Chem. Soc. Faraday Trans. 2 72, 1525-1568.
% 2. Israelachvili, J. N., Marcelja, S. & Horn, R. G. (1980) Q. Rev.
%    Biophys. 13, 121-200.
% 3. Israelachvili, J. N. (1985) Intermolecular and Surface Forces
%    (Academic, New York), pp. 229-271.

% \DLA{This section needs to be added still.}

% PE = 1.33     (Kumar91)
% CHOL = 1.21           (Kumar91)
% PC=0.8        (Kumar91)
% SM=0.8        (assumed by rz same as PC)
% PS=1          (no refs so far; should be close to unity; rz)
% 
% From Epand94LIB: a considerable amount of heat is evolved per mole of
% bilayer stabilizer when such molecules are inserted into membranes
% which are under a high curvature strain. If this energy were coupled
% to a membrane event, such as the conformational change in a protein,
% it could be the driving force responsible for such processes.
% 
% From Epand94LIB: considerable energy may be released upon the
% incorporation of certain molecules into membranes which have a low
% radius of spontaneous curvature. {discuss with DL in terms of
%   catalytic activity; rz; 8/24/2010}.
% 
% Kleinfeld93 showed that the transfer parameters of long-chain FFA
% between the lipid vesicles depend on vesicle curvature and
% composition. Transfer of stearic acid is much slower from LUV as
% compared to SUV.
% 
% Substitution of DPPC SUV for POPC SUV as the donor matrix revealed
% some differences in the energetics of transfer that are associated
% with the physical state of the donor particle. Whereas the free
% energies of transfer of oleic acid from DPPC SUV and POPC SUV are
% similar, the free energy of transfer of stearic acid from DPPC SUV is
% much higher than that observed from POPC.
% 
% From Kumar91: S is given by S = V/al, where V is the hydrocarbon
% volume, a is the area of head group, and l is the critical length of
% the hydrocarbon chain (1-3). a, V, and I are all estimable or
% measurable (4), and the value of S can be calculated. The value of S
% determines the aggregate formed by lipids or any amphiphiles upon
% hydration. It has been shown that lipids aggregate to form spherical
% micelles (S < 1/3), nonspherical (cylindrical) micelles (1/3 < S <
% 1/2), bilayers (1/2 < S < 1), and reverse micelles or hexagonal (HII,)
% phases (S > 1). However, theoreticians caution that the above
% predicted limits set on the values of S are relatively insensitive to
% the exact values of V and a but are strongly dependent upon the choice
% of l (5).
% 
% From Kumar91: Cylindrical and wedge-shaped molecules have been shown
% to be essential for spontaneous vesiculation (6) as well as bilayer
% stabilization (7, 8).
% 
% From Kumar91: It has also been shown experimentally that molecules
% having complementary molecular shapes can form bilayer structures.
% Inverted cone (wedge)- shaped lysoPC and cone-shaped Chol or
% unsaturated PE can interact stoichiometrically to form bilayer
% structures (9-12). {quote to justify our formula for complementary
%   cones; rz}

\section{Kinetic adjustments}
\label{sec:kinetic_adjustments}

In the mass action equation used in the formalism, both the forward
and backwards kinetic constants ($k_\mathrm{f}$ and $k_\mathrm{b}$) are adjusted (by
$k_{\mathrm{f}i\mathrm{adj}}$ and $k_{\mathrm{b}i\mathrm{adj}}$) to reflect the
properties of the vesicle and the properties of the monomer entering
or exiting the vesicle.

\subsection{Forward Kinetic Adjustments}

The forward rate constant adjustment, $k_{\mathrm{f}i\mathrm{adj}}$ takes into
account unsaturation ($un_\mathrm{f}$), charge ($ch_\mathrm{f}$), curvature ($cu_\mathrm{f}$),
length ($l_\mathrm{f}$), and complex formation ($CF1_\mathrm{f}$), each of which is
modified depending on the specific component and the vesicle into
which the component is entering; the final $k_{\mathrm{f}i\mathrm{adj}}$ is the
product of the adjustments made separately for each property.

\begin{equation}
  k_{\mathrm{f}i\mathrm{adj}} = un_\mathrm{f} \cdot ch_\mathrm{f} \cdot cu_\mathrm{f} \cdot l_\mathrm{f} \cdot CF1_\mathrm{f}
  \label{eq:kf_adj}
\end{equation}

\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 heterogeneous 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 %
%%% \RZ{May need to site (at least for us) works showing
%%%   mismatch-dependent ``defects''}. %
Assuming that an increase in width of the distribution linearly
decreases the free energy of activation, the $un_\mathrm{f}$ parameter must
follow $a^{\mathrm{stdev}\left(un_\mathrm{vesicle}\right)}$ where $a > 1$,
so a convenient starting base for $a$ which leads to a reasonable
change in Eyring activation energy is $2$:

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

The mean $\mathrm{stdev}\left(un_\mathrm{vesicle}\right)$ in our
simulations is around $1.5$, which leads to a $\Delta \Delta
G^\ddagger$ of $\Sexpr{to.latex(format(digits=3,to.kcal(2^1.5)))}
\frac{\mathrm{kcal}}{\mathrm{mol}}$, and a total range of possible
values depicted in Figure~\ref{fig:unf_graph}.

% \RZ{Explain here, or even earlier that the formulas were ad hoc
%   adjusted to correspond to ``reasonable'' changes in the Eyring
%   activation energies.}

% 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=6.4in}
\begin{figure}
<<echo=FALSE,results='hide',fig.width=10,fig.height=5>>=
layout(matrix(1:2,nrow=1,ncol=2))
curve(2^x,from=0,to=sd(c(0,4)),
      xlab=expression(paste(stdev,group("(",italic(un[vesicle]),")"))),
      ylab=expression(italic(un[f])))
vps <- baseViewports()
pushViewport(vps$figure)
grid.text("A",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport()
curve(to.kcal(2^x),from=0,to=sd(c(0,4)),
      xlab=expression(paste(stdev,group("(",italic(un[vesicle]),")"))),
      ylab="Activation Energy (kcal/mol)")
vps <- baseViewports()
pushViewport(vps$figure)
grid.text("B",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport()
@ 
\caption{Change in unsaturation forward ($un_\mathrm{f}$) (A) or activation
  energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the
  standard deviation of the unsaturation of the vesicle
  ($\mathrm{stdev}\left(un_\mathrm{vesicle}\right)$).}
\label{fig:unf_graph}
\end{figure}

\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}$ (which is similar to
the numerator of Coulomb's law), 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. The base for $a$ ($60$) was chosen ad hoc to
correspond to a maximum of about
$0.5\frac{\mathrm{kcal}}{\mathrm{mol}}$ change in activation energy
for the common case, resulting in the following equation:


\begin{equation}
  ch_\mathrm{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}}$, and the total range of possible
values seen in Figure~\ref{fig:chf_graph}.


\begin{figure}
<<echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
trellis.device(new=F,color=TRUE)
trellis.par.set(list(axis.line =list(col="transparent")))
pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
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,
                box=NA,
                scales=list(arrows=F,col=1,distance=0.75),
                xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
                ylab=list(label=expression(italic(ch[monomer])),rot=-35),
                zlab=list(label=expression(italic(ch[f])),rot=93)),newpage=FALSE)
rm(x,y,grid)
grid.text("A",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport(1)
pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
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,
                strip=FALSE,
                scales=list(arrows=F,col=1,distance=0.75),
                xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
                ylab=list(label=expression(italic(ch[monomer])),rot=-35),
                zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
rm(x,y,grid)
grid.text("B",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport(2)
@ 
\caption{Change in charge forward ($ch_\mathrm{f}$) (A) or activation energy in
  $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the charge of the
  monomer entering ($ch_\mathrm{monomer}$) and the average charge of
  the vesicle ($\left<ch_\mathrm{vesicle}\right>$).}
\label{fig:chf_graph}
\end{figure}

\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 (neutral
curvature)~\citep{Israelachvili1975:amphiphile_self_assembly}. %
In this formalism, curvature is measured as the ratio of the volume of
the lipid to the area of the head times the length of the lipid
($S=\frac{v}{l_c a}$), so negative curvature is bounded by $(0,1)$,
neutral curvature is 1, and positive curvature is bounded by
$(1,\infty)$. The curvature can be transformed using $\log$, which has
the 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 component with curvature
$cu^{-1}$ will cancel out a component with curvature $cu$, so we have
to log transform (turning these into $-\log cu$ and $\log cu$), then
take the absolute value ($\log cu$ and $\log cu$), and finally measure
the width of the distribution, which in the case of exactly mismatched
curvatures is $0$. 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>}$.
An ad hoc base for $a$ which for the common case is $10$, yielding:


\begin{equation}
 % cu_\mathrm{f} = 10^{\mathrm{stdev}\left|\log cu_\mathrm{vesicle}\right|}
  cu_\mathrm{f} = 10^{\left|\left<\log cu_\mathrm{monomer} \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. The full range of
$cu_\mathrm{f}$ values possible are shown in Figure~\ref{fig:cuf_graph}.

% 1.5 to 0.75 3 to 0.33
\begin{figure}
<<echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
trellis.device(new=F,color=TRUE)
trellis.par.set(list(axis.line =list(col="transparent")))
pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
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,
                box=NA,
                scales=list(arrows=F,col=1,distance=0.75),
                xlab=list(label=expression(paste(stdev,group("(",italic(log(cu[vesicle])),")"))),rot=30),
                ylab=list(label=expression(paste(symbol("\341"),italic(log(cu[vesicle])),symbol("\361"))),rot=-35),
                zlab=list(label=expression(italic(cu[f])),rot=93)),newpage=FALSE)
rm(grid)
grid.text("A",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport(1)
pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
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,
                box=NA,
                scales=list(arrows=F,col=1,distance=0.75),
                xlab=list(label=expression(paste(stdev,group("(",italic(log(cu[vesicle])),")"))),rot=30),
                ylab=list(label=expression(paste(symbol("\341"),italic(log(cu[vesicle])),symbol("\361"))),rot=-35),
                zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
rm(grid)
grid.text("B",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport(2)
@ 
\caption{Change in curvature forward ($cu_\mathrm{f}$) (A) or activation energy
  in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the standard
  deviation of the log of the curvature of the vesicle
  ($\mathrm{stdev}\left(\log cu_\mathrm{vesicle}\right)$) and the mean
  of the log of the curvature of the vesicle ($\left<\log
    cu_\mathrm{vesicle}\right>$).}
\label{fig:cuf_graph}
\end{figure}

\subsubsection{Length Forward}

As in the case of unsaturation, void formation is easier when vesicles
are made up of components of mismatched 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{vesicle}}$, where $\mathrm{stdev}
l_\mathrm{vesicle}$ is the standard deviation of the length of the
components of the vesicle, which has a maximum possible value of
$\Sexpr{format(digits=3,sd(c(rep(12,50),rep(24,50))))}$ and a minimum
of 0 in this set of simulations. Based on activation energy, a
reasonable base for $x$ is 2, leading to:

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

The most common $\mathrm{stdev} l_\mathrm{vesicle}$ 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.

% \DLA{Incorporate McLean84 here}
% From McLean84LIB: Although it is difficult to measure cmc values for
% the sparingly soluble lipids used in this study, estimates for
% lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds,
% 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution;
% Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6
% X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X
% 10-8 M for DMPC was estimated from the linear relationship between ln
% cmc and the number of carbons in the PC acyl chain by using data for n
% = 7, 8, 10, and 16 [summarized in Tanford (1980)].
% 
% From Nichols85: The magnitude of the dissociation rate constant
% decreases by a factor of approximately 3.2 per carbon increase in acyl
% chain length (see Table II here) {take into account for the formula;
%   rz 8/17/2010}.
% 
% From Nichols85: The magnitude of the partition coefficient increases
% with acyl chain length [Keq(M-C6-NBD-PC) = (9.8±2.1) X 106 M-1 and Keq
% (P-C6-NBD-PC) = (9.4±1.3) x 107 M-1. {take into account for the
%   formula; rz 8/17/2010}.
% 
% From Nichols85: The association rate constant is independent of acyl
% chain length. {take into account for the formula; rz 8/17/2010}.



\begin{figure}
<<echo=FALSE,results='hide',fig.width=14,fig.height=5>>=
layout(matrix(1:2,nrow=1,ncol=2))
curve(2^x,from=0,to=sd(c(12,24)),
      xlab=expression(paste(stdev,group("(",italic(l[vesicle]),")"))),
      ylab=expression(italic(l[f])),
      las=1
      )
vps <- baseViewports()
pushViewport(vps$figure)
grid.text("A",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport()
curve(to.kcal(2^x),from=0,to=sd(c(12,24)),
      xlab=expression(paste(stdev,group("(",italic(l[vesicle]),")"))),
      ylab="Activation Energy (kcal/mol)",
      las=1)
vps <- baseViewports()
pushViewport(vps$figure)
grid.text("B",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport()
@ 
\caption{Change in length forward ($l_\mathrm{f}$) (A) or activation energy in
  $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the standard
  deviation of the length of the vesicle
  ($\mathrm{stdev}\left(l_\mathrm{vesicle}\right)$).}
\label{fig:lf_graph}
\end{figure}

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

\begin{equation}
  CF1_\mathrm{f}=1
  \label{eq:complex_formation_forward}
\end{equation}

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

Just as the forward rate constant adjustment $k_{\mathrm{f}i\mathrm{adj}}$
does, the backwards rate constant adjustment $k_{\mathrm{b}i\mathrm{adj}}$
takes into account unsaturation ($un_\mathrm{b}$), charge ($ch_\mathrm{b}$), curvature
($cu_\mathrm{b}$), length ($l_\mathrm{b}$), and complex formation ($CF1_\mathrm{b}$), each of
which are modified depending on the specific component and the vesicle
from which the component is exiting:


\begin{equation}
  k_{\mathrm{b}i\mathrm{adj}} = un_\mathrm{b} \cdot ch_\mathrm{b} \cdot cu_\mathrm{b} \cdot l_\mathrm{b} \cdot CF1_\mathrm{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 (eg. average unsaturation of 2) is more
favorable for lipids with more unsaturation (eg. 3) than lipids with
equivalently less unsaturation (eg. 1), so the difference in energy
between unsaturation is not linear. Therefore, an equation with the
shape $x^{\left| y^{-\left< un_\mathrm{vesicle}\right>
    }-y^{-un_\mathrm{monomer}} \right| }$, where
$\left<un_\mathrm{vesicle}\right>$ is the average unsaturation of the
vesicle and $un_\mathrm{monomer}$ is the unsaturation of the monomer,
allows for increasing the efflux of molecules from membranes where
they strongly mismatch, while allowing vesicles with greater
unsaturation to tolerate greater unsaturation mismatch between
monomers and the vesicle. The bases x and y were chosen ad hoc to
produce reasonable Eyring energies of activation over the range of
unsaturations expected, leading to:


\begin{equation}
  un_\mathrm{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{vesicle}\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. See Figure~\ref{fig:unb_graph} for
the full range of possible values.


\begin{figure}
<<echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
trellis.device(new=F,color=TRUE)
trellis.par.set(list(axis.line =list(col="transparent")))
pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
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,
                box=NA,
                scales=list(arrows=F,col=1,distance=0.75),
                xlab=list(label=expression(paste(symbol("\341"),italic(un[vesicle]),symbol("\361"))),rot=30),
                ylab=list(label=expression(italic(un[monomer])),rot=-35),
                zlab=list(label=expression(italic(un[b])),rot=93)),newpage=FALSE)
grid.text("A",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport(1)
pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
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,
                strip=FALSE,
                scales=list(arrows=F,col=1,distance=0.75),
                xlab=list(label=expression(paste(symbol("\341"),italic(un[vesicle]),symbol("\361"))),rot=30),
                ylab=list(label=expression(italic(un[monomer])),rot=-35),
                zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
rm(grid)
grid.text("B",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport(2)
@ 
\caption{Change in unsaturation backward ($un_\mathrm{b}$) (A) or activation energy in
  $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the unsaturation of the
  monomer leaving ($un_\mathrm{monomer}$) and the average unsaturation of
  the vesicle ($\left<un_\mathrm{vesicle}\right>$).}
\label{fig:unb_graph}
\end{figure}

\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_\mathrm{b} = a^{\left<ch_v\right>
  ch_m}$ is then appropriate, and using a base of $a=20$ yields:

\begin{equation}
  ch_\mathrm{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$.
See Figure~\ref{fig:chb_graph} for the full range of possible values
of $ch_\mathrm{b}$.


\begin{figure}
<<echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
trellis.device(new=F,color=TRUE)
trellis.par.set(list(axis.line =list(col="transparent")))
pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
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,
                box=NA,
                scales=list(arrows=F,col=1,distance=0.75),
                xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
                ylab=list(label=expression(italic(ch[monomer])),rot=-35),
                zlab=list(label=expression(italic(ch[b])),rot=93)),newpage=FALSE)
rm(x,y,grid)
grid.text("A",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport(1)
pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
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,
                strip=FALSE,
                scales=list(arrows=F,col=1,distance=0.75),
                xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
                ylab=list(label=expression(italic(ch[monomer])),rot=-35),
                zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
rm(x,y,grid)
grid.text("B",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport(2)
@ 
\caption{Change in charge backward ($ch_\mathrm{b}$) (A) or activation energy in
  $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the charge of the
  monomer leaving ($ch_\mathrm{monomer}$) and the average charge of
  the vesicle ($\left<ch_\mathrm{vesicle}\right>$).}
\label{fig:chb_graph}
\end{figure}


\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_\mathrm{f}$ needs to be one. To map negative and
positive curvature to the same range, we also need take the logarithm.
Positive and negative curvature lipids cancel each other out,
resulting in an average curvature of 0; the average of the log also
has this property. Increasing mismatches in curvature increase the
rate of efflux, but asymptotically. An equation which satisfies these
criteria has the form $cu_\mathrm{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 absolute value of the difference, however, this
yields a cusp and sharp increase about the curvature equilibrium. We
have chosen bases of $a=7$ and $b=20$ ad hoc, based on activation
energy considerations.

\begin{equation}
  cu_\mathrm{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{vesicle}\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. The full range of values possible for $cu_\mathrm{b}$ are shown in
Figure~\ref{fig:cub_graph}

% \RZ{What about the opposite curvatures that actually do fit to each
%   other?}


% From Kumar91: It has also been shown experimentally that molecules
% having complementary molecular shapes can form bilayer structures.
% Inverted cone (wedge)- shaped lysoPC and cone-shaped Chol or
% unsaturated PE can interact stoichiometrically to form bilayer
% structures (9-12). {quote to justify our formula for complementary
%   cones; rz}


\begin{figure}
<<echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
trellis.device(new=F,color=TRUE)
trellis.par.set(list(axis.line =list(col="transparent")))
pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
grid <- expand.grid(x=log(seq(0.8,1.33,length.out=20)),
                    y=seq(0.8,1.33,length.out=20))
grid$z <- 7^(1-1/(20*(grid$x-log(grid$y))^2+1))
print(wireframe(z~x*y,grid,cuts=50,
                drape=TRUE,
                box=NA,
                scales=list(arrows=F,col=1,distance=0.75),
                xlab=list(label=expression(paste(symbol("\341"),log(italic(cu[vesicle])),symbol("\361"))),rot=30),
                ylab=list(label=expression(italic(cu[monomer])),rot=-35),
                zlab=list(label=expression(italic(cu[b])),rot=93)),newpage=FALSE)
rm(grid)
grid.text("A",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport(1)
pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
grid <- expand.grid(x=log(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*(grid$x-log(grid$y))^2+1)))
print(wireframe(z~x*y,grid,cuts=50,
                drape=TRUE,
                strip=FALSE,
                scales=list(arrows=F,col=1,distance=0.75),
                xlab=list(label=expression(paste(symbol("\341"),log(italic(cu[vesicle])),symbol("\361"))),rot=30),
                ylab=list(label=expression(italic(cu[monomer])),rot=-35),
                zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
rm(grid)
grid.text("B",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport(2)
@ 
\caption{Change in curvature backward ($cu_\mathrm{b}$) (A) or activation
  energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the
  curvature of the monomer leaving ($ch_\mathrm{monomer}$) and the
  average of the log of the curvature of the vesicle
  ($\left<\log cu_\mathrm{vesicle}\right>$).}
\label{fig:cub_graph}
\end{figure}


\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~%
\citep{Nichols1985:phospholipid_monomer_vesicle_thermodynamics}.
Unfortunately, the experimental data known to us only measures the
effect of chain lengths less than or equal to the bulk lipid, not the
effect of lengths exceeding it, and is only known for one bulk lipid
component (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_\mathrm{b} = 3.2^{\left|\left<l_\mathrm{vesicle}\right>-l_\mathrm{monomer}\right|}
  \label{eq:length_backward}
\end{equation}

The most common $\left<l_\mathrm{vesicle}\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 length near 18. The full range of possible values of
$l_\mathrm{b}$ are shown in Figure~\ref{fig:lb_graph}

% (for methods? From McLean84LIB: The activation free energies and free
% energies of transfer from self-micelles to water increase by 2.2 and
% 2.1 kJ mol-' per methylene group, respectively. \RZ{see if we can use it
%   to justify arranging our changed in activating energy around 1
%   kcal/mol}).


\begin{figure}
<<echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
trellis.device(new=F,color=TRUE)
trellis.par.set(list(axis.line =list(col="transparent")))
pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
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))/1e6
print(wireframe(z~x*y,grid,cuts=50,
                drape=TRUE,
                box=NA,
                scales=list(arrows=F,col=1,distance=0.75),
                xlab=list(label=expression(paste(symbol("\341"),italic(l[vesicle]),symbol("\361"))),rot=30),
                ylab=list(label=expression(italic(l[monomer])),rot=-35),
                zlab=list(label=expression(italic(l[b])%*%10^-6),rot=93)),newpage=FALSE)
rm(grid)
grid.text("A",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport(1)
pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
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,
                strip=FALSE,
                scales=list(arrows=F,col=1,distance=0.75),
                xlab=list(label=expression(paste(symbol("\341"),italic(l[vesicle]),symbol("\361"))),rot=30),
                ylab=list(label=expression(italic(l[monomer])),rot=-35),
                zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
rm(grid)
grid.text("B",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport(2)
@ 
\caption{Change in length backwards ($l_\mathrm{b}\times 10^{-6}$) (A) or activation energy
  in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the length of the
  monomer leaving ($l_\mathrm{monomer}$) and the average length of
  lipids in the vesicle ($\left<l_\mathrm{vesicle}\right>$).}
\label{fig:lb_graph}
\end{figure}


\subsubsection{Complex Formation Backward}

Complex formation ($CF1$) describes the interaction between CHOL and
PC or SM, where PC or SM protects the hydroxyl group of CHOL from
interactions with water %
\citep{Huang1999:chol_solubility_model,
  Huang1999:cholesterol_solubility, McConnell2006, McConnell2003}%
. 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), components with negative $CF1$
(CHOL) will be retained. If average $CF1$ is negative, components with
positive $CF1$ are retained. An equation which has this property is
$CF1_\mathrm{b}=a^{\left<CF1_\mathrm{vesicle}\right>
  CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{vesicle}\right>
    CF1_\mathrm{monomer}\right|}$, where difference of the exponent is
zero if the average $CF1$ and the $CF1$ of the component have the same
sign, or double the product if the signs are different. Based on
activation energy considerations, we took an ad hoc base for $a$ as
$1.5$.


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

The most common $\left<CF1_\mathrm{vesicle}\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$. The full range of possible values for $CF1_\mathrm{b}$ are
depicted in Figure~\ref{fig:cf1b_graph}.



\begin{figure}
<<echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
trellis.device(new=F,color=TRUE)
trellis.par.set(list(axis.line =list(col="transparent")))
pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
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,
                box=NA,
                scales=list(arrows=F,col=1,distance=0.75),
                xlab=list(label=expression(paste(symbol("\341"),italic(CF1[vesicle]),symbol("\361"))),rot=30),
                ylab=list(label=expression(italic(CF1[monomer])),rot=-35),
                zlab=list(label=expression(italic(CF1[b])),rot=93)),newpage=FALSE)
rm(grid)
grid.text("A",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport(1)
pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
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,
                strip=FALSE,
                scales=list(arrows=F,col=1,distance=0.75),
                xlab=list(label=expression(paste(symbol("\341"),italic(CF1[vesicle]),symbol("\361"))),rot=30),
                ylab=list(label=expression(italic(CF1[monomer])),rot=-35),
                zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
rm(grid)
grid.text("B",
          just=c("left","top"),
          gp=gpar(fontsize=48),
          y=unit(1,"npc")-unit(0.25,"lines"),
          x=unit(0,"npc")+unit(0.25,"lines"))
popViewport(2)
@ 
\caption{Change in complex formation 1 backwards ($CF1_\mathrm{b}$) (A) or
  activation energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus
  the complex formation of the monomer leaving
  ($CF1_\mathrm{monomer}$) and the average complex formation of the
  vesicle ($\left<CF1_\mathrm{vesicle}\right>$).}
\label{fig:cf1b_graph}
\end{figure}

\section{Simulation Methodology}

\subsection{Overall Architecture}

The simulations are currently run using a program written in perl with
various modules to handle the subsidiary parts. It produces output for
each generation, including the step immediately preceding 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~\citep{R:maincite} (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 (0.71), which is seeded for each run using entropy from
the linux kernel's urandom device.

Code is available upon request.

\subsection{Environment Creation}

\subsubsection{Components}
The environment contains concentrations of components. In the current
set of simulations, there are \Sexpr{1+4*5*7} different possible
components, consisting of PC, PE, PS, SM, and CHOL, with all lipids
except for CHOL having 5 possible unsaturations ranging 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 components
(with the exception of CHOL) represented once until the desired number
of unique components is obtained. CHOL is overrepresented
$\Sexpr{5*7}$ times to be at the level of other lipid types, ensuring
that the probability of CHOL being absent in the environment is the
same as the probability of any one of the other lipid types (PS, PE,
etc.) being absent. This reduces the number of simulations with a
small number of components which are devoid of CHOL.

\subsubsection{Concentrations}
Once the components of the environment have been selected, their
concentrations are 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 at 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 per-simulation 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 probability 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. (For example, if a vesicle was to
have 10 molecules split evenly between three lipid types, the 10th
molecule would be assigned by randomly choosing between the three
lipid types, weighted by their concentration in the environment.)

\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 above until the vesicle splits
forming two progeny vesicles, one of which the program continues to
follow.

\subsubsection{Calculation of Vesicle Properties}

\label{sec:ves_prop}
$S_\mathrm{vesicle}$ 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 assume that the lipid packing is
optimal, and there is no empty space.

For the other vesicle properties (length, unsaturation, charge, and
curvature), we calculate the mean, the standard deviation, the mean of
the log, the mean of the absolute value of the log, the standard
deviation of the log, and the standard deviation of the absolute value
of the log. All cases, when we take the log, we take the log of the
absolute value, and map $\log(0)$ to $0$. For the purpose of plotting
vesicle properties, we only plot the mean of the property.

\subsubsection{Joining and Leaving of Lipid Molecules}

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

$n_i = k_{\mathrm{f}i}k_{\mathrm{f}i_\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]S_\mathrm{vesicle}\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_{\mathrm{b}i}k_{\mathrm{b}i_\mathrm{adj}}C_{i_\mathrm{vesicle}}\mathrm{N_A}dt$.

Before addition or removal of molecules, the properties of the vesicle
are calculated; this is done so that molecules can be considered to be
added and removed at the same instant, even though additions are
handled first programmatically. This also avoids cases where a removal
would have resulted in a negative number of molecules for a particular
type.

\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.

\subsubsection{Vesicle splitting}

If a vesicle has grown to a size which is more than double the
starting vesicle size, the vesicle splits. Molecules are assigned to
the progeny vesicles at random, with each progeny 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 produce later output.

\section{Analyzing output}

The analysis 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}

Two major groups of axes that contribute to vesicle variation between
subsequent generations 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 with different number of molecules to be directly compared
(though it hides any effect of vesicle size).

In order to visualize the variation between and transition of
subsequent generations of vesicles from their initial state in the
simulation, to their final state at the termination of the simulation,
we plot the projection of the generations onto the two principle PCA
axes (and alternatively, any pairing of the axes).

\subsection{Carpet plots}

Carpet plots show the distance/similarity between the composition or
properties of all generations in a single run. The difference between
each group of vesicle can be calculated using Euclidean distance
(properties and compositions) or H similarity (composition only). We
must use Euclidean distance for properties because the H distance
metric is invalid when comparing negative vector elements to positive
vector elements.

In addition to showing distance or similarity, carpet plots also
depict propersomes and composomes as square boxes on the diagonals and
propertypes and compotypes as rectangles off the diagonals, each
propertype or compotype with a distinct color.

\subsection{Propersomes, propertypes, composomes and compotypes}

A generation is considered to be a propersome if it is less than $0.6$
units (by Euclidean distance of normalized properties) away from the
generation immediately following and preceding. Likewise, a generation
is in a composome if its H similarity is more than $0.9$ (by the
normalized H metric) from the preceding and following generations.
Propersomes and composomes are then assigned to propertypes and
compotypes using Paritioning Around Medioids
(PAM). All values of $k$ between 2 and 15
(or the number of propersomes and composomes, if that is less than 15)
are tried, and the clustering with the smallest
silhouette~\citep{Rousseeuw1987:silhouettes} is chosen as the ideal
clustering~\citep{Shenhav2005:pgard}.


\bibliographystyle{unsrtnat}
\bibliography{references.bib}


\end{document}