From a58f8a92449aa803bf4396ae5f62f463ab1caf9c Mon Sep 17 00:00:00 2001 From: don Date: Wed, 27 Jul 2011 23:08:11 +0000 Subject: [PATCH] update notes, start competition work git-svn-id: svn+ssh://hemlock.ucr.edu/srv/svn/misc/trunk/origins_of_life@944 25fa0111-c432-4dab-af88-9f31a2f6ac42 --- Makefile | 2 +- kinetic_formalism.Rnw | 33 +- kinetic_formalism_competition.Rnw | 1579 +++++++++++++++++++++++++++++ references.bib | 1 + to_latex.R | 51 + 5 files changed, 1652 insertions(+), 14 deletions(-) create mode 100644 kinetic_formalism_competition.Rnw create mode 120000 references.bib create mode 100644 to_latex.R diff --git a/Makefile b/Makefile index 7ca3bee..4662459 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,6 @@ #!/bin/make -all: kinetic_formalism.pdf +all: kinetic_formalism.pdf kinetic_formalism_competition.pdf R=R diff --git a/kinetic_formalism.Rnw b/kinetic_formalism.Rnw index b653d81..8c4fe20 100644 --- a/kinetic_formalism.Rnw +++ b/kinetic_formalism.Rnw @@ -62,6 +62,11 @@ <>= 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 @@ -72,19 +77,21 @@ to.kcal <- function(k,temp=300) { \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} - diff --git a/kinetic_formalism_competition.Rnw b/kinetic_formalism_competition.Rnw new file mode 100644 index 0000000..917ef85 --- /dev/null +++ b/kinetic_formalism_competition.Rnw @@ -0,0 +1,1579 @@ +\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 + +<>= +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} + +<>= +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} +<>= +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_m}$ (which is similar to +the numerator of Coulomb's law), where $\left$ is the +average charge of the vesicle, and $ch_m$ is the charge of the +monomer. If either $\left$ 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} +<>= +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$).} +\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} +<>= +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} +<>= +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$ 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} - 2^{-un_\mathrm{monomer}} \right)^2+1\right)^{-1}} + \label{eq:unsaturation_backward} +\end{equation} + +The most common $\left$ 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} +<>= +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$).} +\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_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$ 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} +<>= +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$).} +\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} +<>= +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{monomer}\right|} + \label{eq:length_backward} +\end{equation} + +The most common $\left$ 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} +<>= +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$).} +\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{monomer}-\left|\left + 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{monomer}-\left|\left CF1_\mathrm{monomer}\right|} + \label{eq:complex_formation_backward} +\end{equation} + +The most common $\left$ 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} +<>= +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$).} +\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} diff --git a/references.bib b/references.bib new file mode 120000 index 0000000..434f02c --- /dev/null +++ b/references.bib @@ -0,0 +1 @@ +../paper/references.bib \ No newline at end of file diff --git a/to_latex.R b/to_latex.R new file mode 100644 index 0000000..6858fda --- /dev/null +++ b/to_latex.R @@ -0,0 +1,51 @@ +to.latex <- function(x){ + gsub("\\\\","\\\\\\\\",latexSN(x)) +} +cleanup.tolatex <- function(output) { + gsub("\\\\textrm\\{e\\}(-?)0*(\\d+)","$\\\\times 10^{\\1\\2}$",output); +} + + +table.to.latex <- function(table, + digits=NULL, + format=NULL, + scientific=NULL, + nsmall=0, + colspec="r", + useBooktabs=TRUE, + toprule=if(useBooktabs) "\\toprule" else "\\hline\\hline", + midrule=if(useBooktabs) "\\midrule" else "\\hline", + bottomrule=if(useBooktabs) "\\bottomrule" else "\\hline\\hline",...) { + ans <- "" + header <- c("",colnames(table)); + res <- rbind(header, + cbind(rownames(table), + apply(table,c(1,2),function(x){format(x,digits=digits,nsmall=nsmall,scientific=scientific,...)})#format(table,digits=digits,scientific=scientific) + )) + coefrows <- 2:NROW(res) + res[coefrows,-1] <- sub("(\\*+)","^{\\1}",res[coefrows,-1]) + res[coefrows,-1] <- sub("([eE])\\+?(-?)0*([0-9]+)","$\\\\times 10^{\\2\\3}$", + res[coefrows,-1]) + + tabspec <- rep("l",ncol(res)) + tabspec[-1] <- colspec + tabbegin <- paste("\\begin{tabular}{",paste(tabspec,collapse=""),"}",sep="") + tabend <- "\\end{tabular}" + ans <- c(ans,tabbegin) + if(length(toprule)) + ans <- c(ans,toprule) + for(j in 2:ncol(res)) + ans <- c(ans,paste(" & \\multicolumn{1}{c}{",trimws(res[1,j]),"}",sep="")) + ans <- c(ans,"\\\\") + if(length(midrule)) + ans <- c(ans,midrule) + for(i in 2:NROW(res)) { + ans <- c(ans, + paste(paste(res[i,],collapse=" & "),"\\\\")) + } + if(length(bottomrule)) + ans <- c(ans,bottomrule) + ans <- c(ans,tabend) + structure(ans,class="Latex") + +} -- 2.39.2