--- /dev/null
+\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}
+<<fig=TRUE,echo=FALSE,results=hide,width=10,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}
+<<fig=TRUE,echo=FALSE,results=hide,width=14,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}
+<<fig=TRUE,echo=FALSE,results=hide,width=14,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}
+<<fig=TRUE,echo=FALSE,results=hide,width=14,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}
+<<fig=TRUE,echo=FALSE,results=hide,width=14,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}
+<<fig=TRUE,echo=FALSE,results=hide,width=14,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}
+<<fig=TRUE,echo=FALSE,results=hide,width=14,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}
+<<fig=TRUE,echo=FALSE,results=hide,width=14,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}
+<<fig=TRUE,echo=FALSE,results=hide,width=14,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}
projects / ool / lipid_simulation_formalism.git / commitdiff