\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[hyperfigures,backref,bookmarks,colorlinks]{hyperref} \usepackage[sectionbib,sort&compress,square,numbers]{natbib} \usepackage[margin,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}} \oddsidemargin 0.0in \textwidth 6.5in \raggedbottom \clubpenalty = 10000 \widowpenalty = 10000 \pagestyle{fancy} \author{Don Armstrong} \title{OOL Kinetic Formalisms} %\date{} \onehalfspacing \begin{document} %\maketitle <>= require(lattice) require(grid) # R in cal / mol K to.kcal <- function(k,temp=300) { gasconst <- 1.985 return(-gasconst*temp*log(k)/1000) } @ \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: \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}} \label{eq:state} \end{equation} For $k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]$, $k_{fi}$ has units of $\frac{\mathrm{m}}{\mathrm{s}}$, $k_{fi\mathrm{adj}}$ and $k_{bi\mathrm{adj}}$ are unitless, concentration is in units of $\frac{\mathrm{n}}{\mathrm{L}}$, surface area is in units of $\mathrm{m}^2$, $k_{bi}$ has units of $\frac{1}{\mathrm{s}}$ and $C_{i_\mathrm{ves}}$ has units of $\mathrm{n}$, Thus, we have \begin{equation} \frac{\mathrm{n}}{\mathrm{s}} = \frac{\mathrm{m}}{\mathrm{s}} \frac{\mathrm{n}}{\mathrm{L}} \mathrm{m}^2 \frac{1000\mathrm{L}}{\mathrm{m}^3} - \frac{1}{\mathrm{s}} \mathrm{n} = \frac{\mathrm{m^3}}{\mathrm{s}} \frac{\mathrm{n}}{\mathrm{L}} \frac{1000\mathrm{L}}{\mathrm{m}^3} - \frac{\mathrm{n}}{\mathrm{s}}= \frac{\mathrm{n}}{\mathrm{s}} = 1000 \frac{\mathrm{n}}{\mathrm{s}} - \frac{\mathrm{n}}{\mathrm{s}} \label{eq:state_units} \end{equation} 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 \label{eq:kf_adj} \end{equation} \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} 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. \setkeys{Gin}{width=3.2in} <>= curve(2^x,from=0,to=sd(c(0,4)), main="Unsaturation Forward", xlab="Standard Deviation of Unsaturation of Vesicle", ylab="Unsaturation Forward Adjustment") @ <>= curve(to.kcal(2^x),from=0,to=sd(c(0,4)), main="Unsaturation forward", xlab="Standard Deviation of Unsaturation of Vesicle", ylab="Unsaturation Forward (kcal/mol)") @ \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_m}$, 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. 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} 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}}$. <>= 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), 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) @ <>= 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), 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^{\left|\left<\log cu_\mathrm{vesicle} \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. % 1.5 to 0.75 3 to 0.33 <>= 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) @ <>= 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 = 2^{\mathrm{stdev} l_\mathrm{ves}} \label{eq:length_forward} \end{equation} 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. <>= curve(2^x,from=0,to=sd(c(12,24)), main="Length forward", xlab="Standard Deviation of Length of Vesicle", ylab="Length Forward Adjustment") @ <>= 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} \end{equation} \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: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 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$ 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 = 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. <>= 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, scales=list(arrows=FALSE), xlab=list("Average Vesicle Unsaturation",rot=30), ylab=list("Monomer Unsaturation",rot=-35), zlab=list("Unsaturation Backward",rot=93))) rm(grid) @ <>= 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, scales=list(arrows=FALSE), xlab=list("Average Vesicle Unsaturation",rot=30), ylab=list("Monomer Unsaturation",rot=-35), zlab=list("Unsaturation Backward (kcal/mol)",rot=93))) rm(grid) @ \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_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} 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$. <>= 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, scales=list(arrows=FALSE), xlab=list("Average Vesicle Charge",rot=30), ylab=list("Component Charge",rot=-35), zlab=list("Charge Backwards",rot=93))) rm(x,y,grid) @ <>= 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, scales=list(arrows=FALSE), xlab=list("Average Vesicle Charge",rot=30), ylab=list("Component Charge",rot=-35), zlab=list("Charge Backwards (kcal/mol)",rot=93))) rm(x,y,grid) @ \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_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{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. <>= 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)) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, scales=list(arrows=FALSE), xlab=list("Vesicle Curvature",rot=30), ylab=list("Monomer Curvature",rot=-35), zlab=list("Curvature Backward",rot=93))) rm(grid) @ <>= 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))) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, scales=list(arrows=FALSE), xlab=list("Vesicle Curvature",rot=30), ylab=list("Monomer Curvature",rot=-35), zlab=list("Curvature Backward (kcal/mol)",rot=93))) 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|\left-l_\mathrm{monomer}\right|} \label{eq:length_backward} \end{equation} 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. <>= 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)) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, scales=list(arrows=FALSE), xlab=list("Average Vesicle Length",rot=30), ylab=list("Monomer Length",rot=-35), zlab=list("Length Backward",rot=93))) rm(grid) @ <>= 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, scales=list(arrows=FALSE), xlab=list("Average Vesicle Length",rot=30), ylab=list("Monomer Length",rot=-35), zlab=list("Length Backward (kcal/mol)",rot=93))) rm(grid) @ \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{monomer}-\left|\left 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^{\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$. <>= 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, scales=list(arrows=FALSE), xlab=list("Vesicle Complex Formation",rot=30), ylab=list("Monomer Complex Formation",rot=-35), zlab=list("Complex Formation Backward",rot=93))) rm(grid) @ <>= 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, scales=list(arrows=FALSE), xlab=list("Vesicle Complex Formation",rot=30), ylab=list("Monomer Complex Formation",rot=-35), zlab=list("Complex Formation Backward (kcal/mol)",rot=93))) rm(grid) @ \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} % \bibliography{references.bib} \end{document}