\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
+
+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} -
k_{bi}k_{bi\mathrm{adj}}C_{i_\mathrm{ves}}
The 1000 isn't in \fref{eq:state} above, because it is unit-dependent.
\subsection{Forward adjustments ($k_{fi\mathrm{adj}}$)}
+\label{sec:forward_adj}
+
+The forward rate constant adjustment, $k_{fi\mathrm{adj}}$ takes into
+account unsaturation ($un_f$), charge ($ch_f$), curvature ($cu_f$),
+length ($l_f$), and complex formation ($CF1_f$), each of which are
+modified depending on the specific specie and the vesicle into which
+the specie is entering.
\begin{equation}
k_{fi\mathrm{adj}} = un_f \cdot ch_f \cdot cu_f \cdot l_f \cdot CF1_f
\newpage
\subsubsection{Unsaturation Forward}
+
+In order for a lipid to be inserted into a membrane, a void has to be
+formed for it to fill. Voids can be generated by the combination of
+unsaturated and saturated lipids forming herterogeneous domains. Void
+formation is increased when the unsaturation of lipids in the vesicle
+is widely distributed; in other words, the insertion of lipids into
+the membrane is greater when the standard deviation of the
+unsaturation is larger. Assuming that an increase in width of the
+distribution linearly decreases the free energy of activation, the
+$un_f$ parameter must follow
+$a^{\mathrm{stdev}\left(un_\mathrm{ves}\right)}$ where $a > 1$, so a
+convenient starting base for $a$ is $2$:
+
\begin{equation}
un_f = 2^{\mathrm{stdev}\left(un_\mathrm{ves}\right)}
\label{eq:unsaturation_forward}
\end{equation}
-<<fig=TRUE,echo=FALSE,results=hide,width=5,height=5>>=
+The most common $\mathrm{stdev}\left(un_\mathrm{ves}\right)$ is around
+$1.5$, which leads to a $\Delta \Delta G^\ddagger$ of
+$\Sexpr{format(digits=3,to.kcal(2^1.5))}
+\frac{\mathrm{kcal}}{\mathrm{mol}}$.
+
+It is not clear that the unsaturation of the inserted monomer will
+affect the rate of the insertion positively or negatively, so we do
+not include a term for it in this formalism.
+
+
+<<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",
+ 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",
\newpage
\subsubsection{Charge Forward}
+
+A charged lipid such as PS approaching a vesicle with an average
+charge of the same sign will experience repulsion, whereas those with
+different signs will experience attraction, the degree of which is
+dependent upon the charge of the monomer and the average charge of the
+vesicle. If either the vesicle or the monomer has no charge, there
+should be no effect of charge upon the rate. This leads us to the
+following equation, $a^{-\left<ch_v\right> ch_m}$, where
+$\left<ch_v\right>$ is the average charge of the vesicle, and $ch_m$
+is the charge of the monomer. If either $\left<ch_v\right>$ or $ch_m$
+is 0, the adjustment parameter is 1 (no change), whereas it decreases
+if both are positive or negative, as the product of two real numbers
+with the same sign is always positive. A convenient base for $a$ is
+60, resulting in the following equation:
+
+
\begin{equation}
ch_f = 60^{-\left<{ch}_v\right> {ch}_m}
\label{eq:charge_forward}
\end{equation}
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+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}}$.
+
+<<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)
grid$z <- as.vector(60^(-outer(x,y)))
print(wireframe(z~x*y,grid,cuts=50,
- drape=TRUE,
- scales=list(arrows=FALSE),
- xlab=list("Average Vesicle Charge",rot=30),
- ylab=list("Component Charge",rot=-35),
- zlab=list("Charge Forward",rot=93)))
+ drape=TRUE,
+ scales=list(arrows=FALSE),
+ main="Charge Forward",
+ xlab=list("Average Vesicle Charge",rot=30),
+ ylab=list("Component Charge",rot=-35),
+ zlab=list("Charge Forward",rot=93)))
rm(x,y,grid)
@
-
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+<<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)
grid$z <- as.vector(to.kcal(60^(-outer(x,y))))
print(wireframe(z~x*y,grid,cuts=50,
- drape=TRUE,
- scales=list(arrows=FALSE),
- xlab=list("Average Vesicle Charge",rot=30),
- ylab=list("Component Charge",rot=-35),
- zlab=list("Charge Forward (kcal/mol)",rot=93)))
+ drape=TRUE,
+ scales=list(arrows=FALSE),
+ main="Charge Forward (kcal/mol)",
+ xlab=list("Average Vesicle Charge",rot=30),
+ ylab=list("Component Charge",rot=-35),
+ zlab=list("Charge Forward (kcal/mol)",rot=93)))
rm(x,y,grid)
@
\newpage
\subsubsection{Curvature Forward}
+
+Curvature is a measure of the intrinsic propensity of specific lipids
+to form micelles (positive curvature), inverted micelles (negative
+curvature), or planar sheets (zero curvature). In this formalism,
+curvature is measured as the ratio of the size of the head to that of
+the base, so negative curvature is bounded by $(0,1)$, zero curvature
+is 1, and positive curvature is bounded by $(1,\infty)$. The curvature
+can be transformed into the typical postive/negative mapping using
+$\log$, which has the additional property of making the range of
+positive and negative curvature equal, and distributed about 0.
+
+As in the case of unsaturation, void formation is increased by the
+presence of lipids with mismatched curvature. Thus, a larger
+distribution of curvature in the vesicle increases the rate of lipid
+insertion into the vesicle. However, a species with curvature $e^{-1}$
+will cancel out a species with curvature $e$, so we have to log
+transform (turning these into -1 and 1), then take the absolute value
+(1 and 1), and finally measure the width of the distribution. Thus, by
+using the log transform to make the range of the lipid curvature equal
+between positive and negative, and taking the average to cancel out
+exactly mismatched curvatures, we come to an equation with the shape
+$a^{\left<\log cu_\mathrm{vesicle}\right>}$. A convenient base for $a$
+is $10$, yielding:
+
+
\begin{equation}
- cu_f = 10^{\mathrm{stdev}\left|\log cu_\mathrm{vesicle}\right|}
+ % cu_f = 10^{\mathrm{stdev}\left|\log cu_\mathrm{vesicle}\right|}
+ cu_f = 10^{\left|\left<\log cu_\mathrm{vesicle} \right>\right|\mathrm{stdev} \left|\log cu_\mathrm{vesicle}\right|}
\label{eq:curvature_forward}
\end{equation}
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=5>>=
-curve(10^x,from=0,to=max(c(sd(abs(log(c(0.8,1.33)))),
- sd(abs(log(c(1,1.33)))),
- sd(abs(log(c(0.8,1)))))),
- main="Curvature forward",
- xlab="Standard Deviation of Absolute value of the Log of the Curvature of Vesicle",
- ylab="Curvature Forward Adjustment")
+The most common $\left|\left<\log {cu}_v\right>\right|$ is around
+$0.013$, which with the most common $\mathrm{stdev} \log
+cu_\mathrm{vesicle}$ of $0.213$ leads to a $\Delta \Delta G^\ddagger$
+of $\Sexpr{format(digits=3,to.kcal(10^(0.13*0.213)))}
+\frac{\mathrm{kcal}}{\mathrm{mol}}$. This is a consequence of the
+relatively matched curvatures in our environment.
+
+% 1.5 to 0.75 3 to 0.33
+<<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)),
+ mean(log(c(1,0.33))),
+ mean(log(c(0.33,3)))))),length.out=20))
+grid$z <- 10^(grid$x*grid$y)
+print(wireframe(z~x*y,grid,cuts=50,
+ drape=TRUE,
+ scales=list(arrows=FALSE),
+ xlab=list("Vesicle stdev log curvature",rot=30),
+ ylab=list("Vesicle average log curvature",rot=-35),
+ zlab=list("Vesicle Curvature Forward",rot=93)))
+rm(grid)
@
-
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=5>>=
-curve(to.kcal(10^x),from=0,to=max(c(sd(abs(log(c(0.8,1.33)))),
- sd(abs(log(c(1,1.33)))),
- sd(abs(log(c(0.8,1)))))),
- main="Curvature forward",
- xlab="Standard Deviation of Absolute value of the Log of the Curvature of Vesicle",
- ylab="Curvature Forward Adjustment (kcal/mol)")
+<<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)),
+ mean(log(c(1,0.33))),
+ mean(log(c(0.33,3)))))),length.out=20))
+grid$z <- to.kcal(10^(grid$x*grid$y))
+print(wireframe(z~x*y,grid,cuts=50,
+ drape=TRUE,
+ scales=list(arrows=FALSE),
+ xlab=list("Vesicle stdev log curvature",rot=30),
+ ylab=list("Vesicle average log curvature",rot=-35),
+ zlab=list("Vesicle Curvature Forward (kcal/mol)",rot=93)))
+rm(grid)
@
-
\newpage
\subsubsection{Length Forward}
+
+As in the case of unsaturation, void formation is easier when vesicles
+are made up of components of widely different lengths. Thus, when the
+width of the distribution of lengths is larger, the forward rate
+should be greater as well, leading us to an equation of the form
+$x^{\mathrm{stdev} l_\mathrm{ves}}$, where $\mathrm{stdev}
+l_\mathrm{ves}$ is the standard deviation of the length of the
+components of the vesicle, which has a maximum possible value of 8 and
+a minimum of 0 in this set of experiments. A convenient base for $x$
+is 2, leading to:
+
\begin{equation}
- l_f = 3^{\mathrm{stdev} l_\mathrm{ves}}
+ l_f = 2^{\mathrm{stdev} l_\mathrm{ves}}
\label{eq:length_forward}
\end{equation}
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=5>>=
-curve(3^x,from=0,to=sd(c(12,24)),
+The most common $\mathrm{stdev} l_\mathrm{ves}$ is around $3.4$, which leads to
+a range of $\Delta \Delta G^\ddagger$ of
+$\Sexpr{format(digits=3,to.kcal(2^(3.4)))}
+\frac{\mathrm{kcal}}{\mathrm{mol}}$.
+
+While it could be argued that increased length of the monomer could
+affect the rate of insertion into the membrane, it's not clear whether
+it would increase (by decreasing the number of available hydrogen
+bonds, for example) or decrease (by increasing the time taken to fully
+insert the acyl chain, for example) the rate of insertion or to what
+degree, so we do not take it into account in this formalism.
+
+\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}.
+
+
+<<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>>=
-curve(to.kcal(3^x),from=0,to=sd(c(12,24)),
+<<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",
ylab="Length Forward Adjustment (kcal/mol)")
\subsubsection{Complex Formation}
+There is no contribution of complex formation to the forward reaction
+rate in the current formalism.
+
\begin{equation}
CF1_f=1
\label{eq:complex_formation_forward}
\subsection{Backward adjustments ($k_{bi\mathrm{adj}}$)}
+Just as the forward rate constant adjustment $k_{fi\mathrm{adj}}$
+does, the backwards rate constant adjustment $k_{bi\mathrm{adj}}$
+takes into account unsaturation ($un_b$), charge ($ch_b$), curvature
+($cu_b$), length ($l_b$), and complex formation ($CF1_b$), each of
+which are modified depending on the specific specie and the vesicle
+into which the specie is entering:
+
+
\begin{equation}
k_{bi\mathrm{adj}} = un_b \cdot ch_b \cdot cu_b \cdot l_b \cdot CF1_b
- \label{eq:kf_adj}
+ \label{eq:kb_adj}
\end{equation}
-\newpage
\subsubsection{Unsaturation Backward}
+
+Unsaturation also influences the ability of a lipid molecule to leave
+a membrane. If a molecule has an unsaturation level which is different
+from the surrounding membrane, it will be more likely to leave the
+membrane. The more different the unsaturation level is, the greater
+the propensity for the lipid molecule to leave. However, a vesicle
+with some unsaturation is more favorable for lipids with more
+unsaturation than the equivalent amount of less unsatuturation, so the
+difference in energy between unsaturation is not linear. Therefore, an
+equation with the shape
+$x^{\left| y^{-\left< un_\mathrm{ves}\right> }-y^{-un_\mathrm{monomer}} \right| }$
+where $\left<un_\mathrm{ves}\right>$ is the average unsaturation of
+the vesicle, and $un_\mathrm{monomer}$ is the average unsaturation. In
+this equation, as the average unsaturation of the vesicle is larger,
+
\begin{equation}
- un_b = 10^{\left|3.5^{-\left<un_\mathrm{ves}\right>}-3.5^{-un_\mathrm{monomer}}\right|}
+ un_b = 7^{1-\left(20\left(2^{-\left<un_\mathrm{vesicle} \right>} - 2^{-un_\mathrm{monomer}} \right)^2+1\right)^{-1}}
\label{eq:unsaturation_backward}
\end{equation}
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+The most common $\left<un_\mathrm{ves}\right>$ is around $1.7$, which leads to
+a range of $\Delta \Delta G^\ddagger$ from
+$\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-0)^2+1))))}
+\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with 0 unsaturation
+to
+$\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-4)^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
+for monomers with 4 unsaturations.
+
+
+<<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 <- 10^(abs(3.5^-grid$x-3.5^-grid$y))
+grid$z <- (7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1)))
print(wireframe(z~x*y,grid,cuts=50,
drape=TRUE,
scales=list(arrows=FALSE),
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(10^(abs(3.5^-grid$x-3.5^-grid$y)))
+grid$z <- to.kcal((7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1))))
print(wireframe(z~x*y,grid,cuts=50,
drape=TRUE,
scales=list(arrows=FALSE),
@
+
\newpage
\subsubsection{Charge Backwards}
+As in the case of monomers entering a vesicle, monomers leaving a
+vesicle leave faster if their charge has the same sign as the average
+charge vesicle. An equation of the form $ch_b = a^{\left<ch_v\right>
+ ch_m}$ is then appropriate, and using a base of $a=20$ yields:
+
\begin{equation}
ch_b = 20^{\left<{ch}_v\right> {ch}_m}
\label{eq:charge_backwards}
\end{equation}
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+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$.
+
+
+<<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)
\newpage
\subsubsection{Curvature Backwards}
+
+The less a monomer's intrinsic curvature matches the average curvature
+of the vesicle in which it is in, the greater its rate of efflux. If
+the difference is 0, $cu_f$ needs to be one. To map negative and
+positive curvature to the same range, we also need take the logarithm.
+Increasing mismatches in curvature increase the rate of efflux, but
+asymptotically. An equation which satisfies this critera has the
+form $cu_f = a^{1-\left(b\left( \left< \log cu_\mathrm{vesicle} \right>
+ -\log cu_\mathrm{monomer}\right)^2+1\right)^{-1}}$. An
+alternative form would use the aboslute value of the difference,
+however, this yields a cusp and sharp increase about the curvature
+equilibrium, which is decidedly non-elegant. We have chosen bases of
+$a=7$ and $b=20$.
+
\begin{equation}
- cu_f = 7^{1-\left(20\left(\log_{e} cu_\mathrm{vesicle}-\log_{e} cu_\mathrm{monomer}\right)^2+1\right)^{-1}}
+ cu_b = 7^{1-\left(20\left(\left< \log cu_\mathrm{vesicle} \right> -\log cu_\mathrm{monomer} \right)^2+1\right)^{-1}}
\label{eq:curvature_backwards}
\end{equation}
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+The most common $\left<\log cu_\mathrm{ves}\right>$ is around $-0.013$, which leads to
+a range of $\Delta \Delta G^\ddagger$ from
+$\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(0.8))^2+1))))}
+\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature 0.8
+to
+$\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(1.3))^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
+for monomers with curvature 1.3 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature near 1.
+
+
+<<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)))
rm(grid)
@
-
\newpage
\subsubsection{Length Backwards}
+
+In a model membrane, the dissociation constant increases by a factor
+of approximately 3.2 per carbon decrease in acyl chain length (Nichols
+1985). Unfortunatly, the experimental data known to us only measures
+chain length less than or equal to the bulk lipid, and does not exceed
+it, and is only known for one bulk lipid species (DOPC). We assume
+then, that the increase is in relationship to the average vesicle, and
+that lipids with larger acyl chain length will also show an increase
+in their dissociation constant.
+
\begin{equation}
- l_b = 3.2^{\left|l_\mathrm{ves}-l_\mathrm{monomer}\right|}
+ l_b = 3.2^{\left|\left<l_\mathrm{ves}\right>-l_\mathrm{monomer}\right|}
\label{eq:length_backward}
\end{equation}
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+The most common $\left<\log l_\mathrm{ves}\right>$ is around $17.75$, which leads to
+a range of $\Delta \Delta G^\ddagger$ from
+$\Sexpr{format(digits=3,to.kcal(3.2^abs(12-17.75)))}
+\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with length 12
+to
+$\Sexpr{format(digits=3,to.kcal(3.2^abs(24-17.75)))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
+for monomers with length 24 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature near 18.
+
+
+<<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)))
\newpage
\subsubsection{Complex Formation Backward}
+
+Complex formation describes the interaction between CHOL and PC or SM,
+where PC or SM protects the hydroxyl group of CHOL from interactions
+with water, the ``Umbrella Model''. PC ($CF1=2$) can interact with two
+CHOL, and SM ($CF1=3$) with three CHOL ($CF1=-1$). If the average of
+$CF1$ is positive (excess of SM and PC with regards to complex
+formation), species with negative $CF1$ (CHOL) will be retained. If
+average $CF1$ is negative, species with positive $CF1$ are retained.
+An equation which has this property is
+$CF1_b=a^{\left<CF1_\mathrm{ves}\right>
+ CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{ves}\right>
+ CF1_\mathrm{monomer}\right|}$, where difference of the exponent is
+zero if the average $CF1$ and the $CF1$ of the specie have the same
+sign, or double the product if the signs are different. A convenient
+base for $a$ is $1.5$.
+
+
\begin{equation}
- CF1_b=1.5^{CF1_\mathrm{ves} CF1_\mathrm{monomer}-\left|CF1_\mathrm{ves} CF1_\mathrm{monomer}\right|}
+ CF1_b=1.5^{\left<CF1_\mathrm{ves}\right> CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{ves}\right> CF1_\mathrm{monomer}\right|}
\label{eq:complex_formation_backward}
\end{equation}
-<<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
+The most common $\left<CF1_\mathrm{ves}\right>$ is around $0.925$,
+which leads to a range of $\Delta \Delta G^\ddagger$ from
+$\Sexpr{format(digits=3,to.kcal(1.5^(0.925*-1-abs(0.925*-1))))}
+\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
+formation $-1$ to
+$\Sexpr{format(digits=3,to.kcal(1.5^(0.925*2-abs(0.925*2))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
+for monomers with complex formation $2$ to
+$0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
+formation $0$.
+
+
+<<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 <- 3.2^(grid$x*grid$y-abs(grid$x*grid$y))
+grid$z <- 1.5^(grid$x*grid$y-abs(grid$x*grid$y))
print(wireframe(z~x*y,grid,cuts=50,
drape=TRUE,
scales=list(arrows=FALSE),
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(3.2^(grid$x*grid$y-abs(grid$x*grid$y)))
+grid$z <- to.kcal(1.5^(grid$x*grid$y-abs(grid$x*grid$y)))
print(wireframe(z~x*y,grid,cuts=50,
drape=TRUE,
scales=list(arrows=FALSE),
+\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}
+
+The simulation is currently run by single program written in perl
+using various different modules to handle the subsidiary parts. It
+produces output for each generation, including the step immediately
+preceeding and immediately following a vesicle split (and optionally,
+each step) that is written to a state file which contains the state of
+the vesicle, environment, kinetic parameters, program invocation
+options, time, and various other parameters necessary to recreate the
+state vector at a given time. This output file is then read by a
+separate program in perl to produce different output which is
+generated out-of-band; output which includes graphs and statistical
+analysis is performed using R (and various grid graphics modules)
+called from the perl program.
+
+The separation of simulation and output generation allows refining
+output, and simulations performed on different versions of the
+underlying code to be compared using the same analysis methods and
+code. It also allows later simulations to be restarted from a specific
+generation, utilizing the same environment.
+
+Random variables of different distributions are calculated using
+Math::Random, which is seeded for each run using entropy from the
+linux kernel's urandom device.
+
+\subsection{Environment Creation}
+
+\subsubsection{Components}
+The environment contains different concentrations of different
+components. In the current set of simulations, there are
+\Sexpr{1+4*5*7} different components, consisting of PC, PE, PS, SM,
+and CHOL, with all lipids except for CHOL having 5 possible
+unsaturations rangiong from 0 to 4, and 7 lengths from $12,14,...,22$
+($7\cdot 5\cdot4+1=\Sexpr{1+4*5*7}$). In cases where the environment
+has less than the maximum number of components, the components are
+selected in order without replacement from a randomly shuffled deck of
+component (with the exception of CHOL) represented once until the
+desired number of unique components are obtained. CHOL is over
+representated $\Sexpr{5*7}$ times to be at the level of other lipid
+types, ensures that the probability of CHOL being asbent in the
+environment is the same as the probability of one of the other lipid
+types (PS, PE, etc.) being entirely absent. This reduces the number of
+simulations with a small number of components which are entirely
+devoid of CHOL.
+
+\subsubsection{Concentration}
+Once the components of the environment have been selected, the
+concentration of thoes components is determined. In experiments where
+the environmental concentration is the same across all lipid
+components, the concentration is set to $10^{-10}\mathrm{M}$. In other
+cases, the environmental concentration is set to a random number from
+a gamma distribution with shape parameter 2 and an average of
+$10^{-10}\mathrm{M}$. The base concentration ($10^{-10}\mathrm{M}$)
+can also be altered in the initialization of the experiment to
+specific values for specific lipid types or components.
+
+\subsection{Initial Vesicle Creation}
+
+Initial vesicles are seeded by selecting lipid molecules from the
+environment until the vesicle reaches a specific starting size. The
+vesicle starting size has gamma distribution with shape parameter 2
+and a mean of the experimentally specified starting size, with a
+minimum of 5 lipid molecules. Lipid molecules are then selected to be
+added to the lipid membrane according to four different methods. In
+the constant method, molecules are added in direct proportion to their
+concentration in the environment. The uniform method adds molecules in
+proportion to their concentration in the environment scaled by a
+uniform random value, whereas the random method adds molecules in
+proportion to a uniform random value. The final method is a binomial
+method, which adjust the porbability of adding a molecule of a
+specific component by the concentration of that component, and then
+adds the molecules one by one to the membrane. This final method is
+also used in the first three methods to add any missing molecules to
+the starting vesicle which were unallocated due to the requirement to
+add integer numbers of molecules.
+
+\subsection{Simulation Step}
+
+Once the environment has been created and the initial vesicle has been
+formed, molecules join and leave the vesicle based on the kinetic
+parameters and state equation discussed until the vesicle splits
+forming two daughter vesicles, one of which the program continues to
+follow.
+
+\subsubsection{Calculation of Vesicle Properties}
+
+\label{sec:ves_prop}
+$S_\mathrm{ves}$ is the surface area of the vesicle, and is the sum of
+the surface area of all of the individual lipid components; each lipid
+type has a different surface area; we na\"ively assume that the lipid
+packing is optimal, and there is no empty space.
+
+\subsubsection{Joining and Leaving of Lipid Molecules}
+
+Determining the number of molecules to add to the lipid membrane
+($n_i$) requires knowing $k_{fi\mathrm{adj}}$, the surface area of the
+vesicle $S_\mathrm{ves}$ (see~\fref{sec:ves_prop}), the time interval
+$dt$ during which lipids are added, the base $k_{fi}$, and the
+concentration of the monomer in the environment
+$\left[C_{i\mathrm{ves}}\right]$ (see~\fref{eq:state}).
+$k_{fi\mathrm{adj}}$ is calculated (see~\fref{eq:kf_adj}) based on the
+vesicle properties and their hypothesized effect on the rate (in as
+many cases as possible, experimentally based)
+(see~\fref{sec:forward_adj} for details). $dt$ can be varied
+(see~\fref{sec:step_duration}), but for a given step is constant. This
+leads to the following:
+
+$n_i = k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]S_\mathrm{ves}\mathrm{N_A}dt$
+
+In the cases where $n_i > 1$, the integer number of molecules is
+added. Fractional $n_i$ or the fractional remainder after the addition
+of the integer molecules are the probability of adding a specific
+molecule, and are compared to a uniformly distributed random value
+between 0 and 1. If the random value is less than or equal to the
+fractional part of $n_i$, an additional molecule is added.
+
+Molecules leaving the vesicle are handled in a similar manner, with
+
+$n_i = k_{bi}k_{bi\mathrm{adj}}C_{i_\mathrm{ves}}\mathrm{N_A}dt$.
+
+While programatically, the molecule removal happens after the
+addition, the properties that each operates on are the same, so they
+can be considered to have been added and removed at the same instant.
+[This also avoids cases where a removal would have resulted in
+negative molecules for a particular type, which is obviously
+disallowed.]
+
+\subsubsection{Step duration}
+\label{sec:step_duration}
+
+$dt$ is the time taken for each step of joining, leaving, and checking
+split. In the current implementation, it starts with a value of
+$10^{-6}\mathrm{s}$ but this is modified in between each step if the
+number of molecules joining or leaving is too large or small. If more
+than half of the starting vesicle size molecules join or leave in a
+single step, $dt$ is reduced by half. If less than the starting
+vesicle size molecules join or leave in 100 steps, $dt$ is doubled.
+This is necessary to curtail run times and to automatically adjust to
+different experimental runs.
+
+(In every run seen so far, the initial $dt$ was too small, and was
+increased before the first generation occured; at no time was $dt$ too
+large.)
+
+\subsubsection{Vesicle splitting}
+
+If a vesicle has grown to a size which is more than double the
+starting vesicle size, the vesicle splits. More elaborate mechanisms
+for determining whether a vesicle should split are of course possible,
+but not currently implemented. Molecules are assigned to the daughter
+vesicles at random, with each daughter vesicle having an equal chance
+of getting a single molecule. The number of molecules to assign to the
+first vesicle has binomial distribution with a probability of an event
+in each trial of 0.5 and a number of trials equal to the number of
+molecules.
+
+\subsection{Output}
+
+The environment, initial vesicle, and the state of the vesicle
+immediately before and immediately after splitting are stored to disk
+to produce later output.
+
+\section{Analyzing output}
+
+Analyzing of output is handled by a separate perl program which shares
+many common modules with the simulation program. Current output
+includes simulation progress, summary tables, summary statistics, and
+various graphs.
+
+\subsection{PCA plots}
+
+Vesicles have many different axes which contribute to their variation
+between subsequent generations; two major groups of axes are the
+components and properties of vesicles. Each component in a vesicle is
+an axis on its own; it can be measured either as an absolute number of
+molecules in each component, or the fraction of molecules of that
+component over the total number of molecules; the second approach is
+often more convenient, as it allows vesicles of different number of
+molecules to be more directly compared (though it hides any affect of
+vesicle size).
+
+In order to visualize the transition of subsequent generations of
+vesicles from their initial state in the simulation, to their final
+state at the termination of
+
+\subsection{Carpet plots}
+
% \bibliographystyle{plainnat}
projects / ool / lipid_simulation_formalism.git / blobdiff