\begin{document}
%\maketitle
-<<results=hide,echo=FALSE>>=
-require(lattice)
-require(grid)
+<<load_libraries,results="hide",message=FALSE,warning=FALSE,echo=FALSE>>=
+require("lattice")
+require("grid")
+require("Hmisc")
+require("gridBase")
+opts_chunk$set(dev="cairo_pdf",
+ out.width="0.8\\textwidth",
+ out.height="0.8\\textheight",
+ out.extra="keepaspectratio")
+opts_chunk$set(cache=TRUE, autodep=TRUE)
+options(device = function(file, width = 6, height = 6, ...) {
+ cairo_pdf(tempfile(), width = width, height = height, ...)
+ })
+to.latex <- function(x){
+ gsub("\\\\","\\\\\\\\",latexSN(x))
+}
# R in cal / mol K
to.kcal <- function(k,temp=300) {
gasconst <- 1.985
\section{State Equation}
% double check this with the bits in the paper
-Given a base forward kinetic parameter for the $i$th specie $k_{fi}$
-(which is dependent on lipid type, that is PC, PE, PS, etc.), an
-adjustment parameter $k_{fi\mathrm{adj}}$ based on the vesicle and the
-specific specie (length, unsaturation, etc.) (see~\fref{eq:kf_adj}),
-the molar concentration of monomer of the $i$th specie
-$\left[C_{i_\mathrm{monomer}}\right]$, the surface area of the vesicle
-$S_\mathrm{ves}$, the base backwards kinetic parameter for the $i$th
-specie $k_{bi}$ which is also dependent on lipid type, its adjustment
-parameter $k_{bi\mathrm{adj}}$ (see~\fref{eq:kb_adj}), and the molar
-concentration of the $i$th specie in the vesicle
-$\left[C_{i_\mathrm{ves}}\right]$, the change in concentration of the
-$i$th specie in the vesicle per change in time $\frac{d
- C_{i_\mathrm{ves}}}{dt}$ can be calculated:
+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.
\begin{equation}
\frac{d C_{i_\mathrm{ves}}}{dt} = k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]S_\mathrm{ves} -
not include a term for it in this formalism.
-\setkeys{Gin}{width=3.2in}
-<<fig=TRUE,echo=FALSE,results=hide,width=5,height=5>>=
+<<echo=FALSE,results="hide",fig.width=5,fig.height=5,out.width="3.2in">>=
curve(2^x,from=0,to=sd(c(0,4)),
main="Unsaturation Forward",
xlab="Standard Deviation of Unsaturation of Vesicle",
ylab="Unsaturation Forward Adjustment")
@
-<<fig=TRUE,echo=FALSE,results=hide,width=5,height=5>>=
+<<echo=FALSE,results="hide",fig.width=5,fig.height=5,out.width="3.2in">>=
curve(to.kcal(2^x),from=0,to=sd(c(0,4)),
main="Unsaturation forward",
xlab="Standard Deviation of Unsaturation of Vesicle",
$\Sexpr{format(digits=3,to.kcal(60^(-.165*-1)))}
\frac{\mathrm{kcal}}{\mathrm{mol}}$ to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$.
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=7>>=
x <- seq(-1,0,length.out=20)
y <- seq(-1,0,length.out=20)
grid <- expand.grid(x=x,y=y)
zlab=list("Charge Forward",rot=93)))
rm(x,y,grid)
@
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=7>>=
x <- seq(-1,0,length.out=20)
y <- seq(-1,0,length.out=20)
grid <- expand.grid(x=x,y=y)
relatively matched curvatures in our environment.
% 1.5 to 0.75 3 to 0.33
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=7>>=
grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
sd(abs(log(c(1,0.33)))),sd(abs(log(c(0.33,3)))))),length.out=20),
y=seq(0,max(c(mean(log(c(1,3)),
zlab=list("Vesicle Curvature Forward",rot=93)))
rm(grid)
@
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=7>>=
grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
sd(abs(log(c(1,0.33)))),sd(abs(log(c(0.33,3)))))),length.out=20),
y=seq(0,max(c(mean(log(c(1,3)),
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.
+\fixme{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}.
+
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=5>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=5>>=
curve(2^x,from=0,to=sd(c(12,24)),
main="Length forward",
xlab="Standard Deviation of Length of Vesicle",
ylab="Length Forward Adjustment")
@
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=5>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=5>>=
curve(to.kcal(2^x),from=0,to=sd(c(12,24)),
main="Length forward",
xlab="Standard Deviation of Length of Vesicle",
for monomers with 4 unsaturations.
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=7>>=
grid <- expand.grid(x=seq(0,4,length.out=20),
y=seq(0,4,length.out=20))
grid$z <- (7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1)))
zlab=list("Unsaturation Backward",rot=93)))
rm(grid)
@
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=7>>=
grid <- expand.grid(x=seq(0,4,length.out=20),
y=seq(0,4,length.out=20))
grid$z <- to.kcal((7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1))))
for monomers with charge $0$.
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=7>>=
x <- seq(-1,0,length.out=20)
y <- seq(-1,0,length.out=20)
grid <- expand.grid(x=x,y=y)
zlab=list("Charge Backwards",rot=93)))
rm(x,y,grid)
@
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=7>>=
x <- seq(-1,0,length.out=20)
y <- seq(-1,0,length.out=20)
grid <- expand.grid(x=x,y=y)
for monomers with curvature 1.3 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature near 1.
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=7>>=
grid <- expand.grid(x=seq(0.8,1.33,length.out=20),
y=seq(0.8,1.33,length.out=20))
grid$z <- 7^(1-1/(20*(log(grid$x)-log(grid$y))^2+1))
zlab=list("Curvature Backward",rot=93)))
rm(grid)
@
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=7>>=
grid <- expand.grid(x=seq(0.8,1.33,length.out=20),
y=seq(0.8,1.33,length.out=20))
grid$z <- to.kcal(7^(1-1/(20*(log(grid$x)-log(grid$y))^2+1)))
for monomers with length 24 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature near 18.
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=7>>=
grid <- expand.grid(x=seq(12,24,length.out=20),
y=seq(12,24,length.out=20))
grid$z <- 3.2^(abs(grid$x-grid$y))
zlab=list("Length Backward",rot=93)))
rm(grid)
@
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=7>>=
grid <- expand.grid(x=seq(12,24,length.out=20),
y=seq(12,24,length.out=20))
grid$z <- to.kcal(3.2^(abs(grid$x-grid$y)))
formation $0$.
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=7>>=
grid <- expand.grid(x=seq(-1,3,length.out=20),
y=seq(-1,3,length.out=20))
grid$z <- 1.5^(grid$x*grid$y-abs(grid$x*grid$y))
zlab=list("Complex Formation Backward",rot=93)))
rm(grid)
@
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+<<echo=FALSE,results="hide",fig.width=7,fig.height=7>>=
grid <- expand.grid(x=seq(-1,3,length.out=20),
y=seq(-1,3,length.out=20))
grid$z <- to.kcal(1.5^(grid$x*grid$y-abs(grid$x*grid$y)))
@
+
+\subsection{Per-Lipid Kinetic Parameters}
+
+Each of the 5 lipid types have different kinetic parameters; to the
+greatest extent possible, we have derived these from literature.
+
+\begin{table}
+ \centering
+ \begin{tabular}{c c c c c c c}
+ Type & $k_f$ & $k_b$ & Area (\r{A}$^2$) & Charge & CF1 & Curvature \\
+ \hline
+ PC & $3.7\cdot 10^6$ & $2\cdot 10^{-5}$ & 63 & 0 & 2 & 0.8 \\
+ PS & $3.7\cdot 10^6$ & $1.5\cdot 10^{-5}$ & 54 & -1 & 0 & 1 \\
+ CHOL & $5.1\cdot 10^7$ & $2.8\cdot 10^{-4}$ & 38 & 0 & -1 & 1.21 \\
+ SM & $3.7\cdot 10^6$ & $3.1\cdot 10^{-3}$ & 51 & 0 & 3 & 0.8 \\
+ PE & $2.3\cdot 10^6$ & $10^{-5}$ & 55 & 0 & 0 & 1.33 \\
+ \end{tabular}
+ \caption{Kinetic parameters of lipid types}
+ \label{tab:kinetic_parameters_lipid_types}
+\end{table}
+
+\subsubsection{$k_f$ for lipid types}
+For PC, $k_f$ was measured by Nichols85 to be $3.7\cdot 10^6
+\frac{1}{\mathrm{M}\cdot \mathrm{s}}$ by the partitioning of
+P-C$_6$-NBD-PC between DOPC vesicles and water. The method utilized by
+Nichols85 has the weakness of using NBD-PC, with associated label
+perturbations. As similar measures do not exist for SM or PS, we
+assume that they have the same $k_f$. For CHOL, Estronca07 found a
+value for $k_f$ of $5.1\cdot 10^7 \frac{1}{\mathrm{M}\cdot
+ \mathrm{s}}$. For PE, Abreu04 found a value for $k_f$ of $2.3\cdot
+10^6$. \fixme{I'm missing the notes on these last two papers, so this
+ isn't correct yet.}
+
+\subsubsection{$k_b$ for lipid types}
+
+$k_b$ for PC was measured by Wimley90 using a radioactive label and
+large unilammelar vesicles at 30\textdegree C. The other values were
+calculated from the experiments of Nichols82 where the ratio of $k_b$
+of different types was measured to that of PC.
+See~\fref{tab:kinetic_parameters_lipid_types}.
+
+assigned accordingly. kb(PS) was assumed to be the same as kb(PG)
+given by Nichols82 (also ratio from kb(PC)). kb(SM) is taken from
+kb(PC) of Wimley90 (radioactive), and then a ratio of kb(PC)/kb(SM)
+taken from Bai97: = 34/2.2 = 15.45; 2.0 x 10-4 x 15.45 = 3.1 x 10-3 s
+-1. kb(CHOL) taken from Jones90 (radioactive; POPC LUV; 37°).
+
+PC 0.89
+PE 0.45 <- from Nichols82
+PG=PS
+
+
+kb PC is from table 2 of Wimley90, where we have a half life of 9.6
+hours for DMPC. \Sexpr{log(2)/(9.6*60*60)}.
+
+
+
+\subsubsection{Area for lipid types}
+
+
+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}
+
+From Thomas88a: SM decreases the rate of cholesterol transfer, while
+phosphatidylethanolamine (PE) and phosphatidylserine (PS) have no
+effect at physiologically significant levels.
+
+
\section{Simulation Methodology}
\subsection{Overall Architecture}
(see~\fref{sec:step_duration}), but for a given step is constant. This
leads to the following:
-$n_i = k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]S_\mathrm{ves}dt\mathrm{NA}$
+$n_i = k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]S_\mathrm{ves}\mathrm{N_A}dt$
In the cases where $n_i > 1$, the integer number of molecules is
added. Fractional $n_i$ or the fractional remainder after the addition
Molecules leaving the vesicle are handled in a similar manner, with
-$n_i = k_{bi}k_{bi\mathrm{adj}}C_{i_\mathrm{ves}}dt\mathrm{NA}$.
+$n_i = k_{bi}k_{bi\mathrm{adj}}C_{i_\mathrm{ves}}\mathrm{N_A}dt$.
While programatically, the molecule removal happens after the
addition, the properties that each operates on are the same, so they
\section{Analyzing output}
+Analyzing of output is handled by a separate perl program which shares
+many common modules with the simulation program. Current output
+includes simulation progress, summary tables, summary statistics, and
+various graphs.
+
\subsection{PCA plots}
+Vesicles have many different axes which contribute to their variation
+between subsequent generations; two major groups of axes are the
+components and properties of vesicles. Each component in a vesicle is
+an axis on its own; it can be measured either as an absolute number of
+molecules in each component, or the fraction of molecules of that
+component over the total number of molecules; the second approach is
+often more convenient, as it allows vesicles of different number of
+molecules to be more directly compared (though it hides any affect of
+vesicle size).
+
+In order to visualize the transition of subsequent generations of
+vesicles from their initial state in the simulation, to their final
+state at the termination of
+
\subsection{Carpet plots}