X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=kinetic_formalism.Rnw;h=fccc4af7c3a6eb4be8a85ab899b553b8453d4e8c;hb=4ac6dcb5960639dce3d89699983f55a90e022636;hp=8b387e28d27c3ae07b3de4e4235f084dab454267;hpb=b09b39a7a1bb929470680062a8ccbab0da67590d;p=ool%2Flipid_simulation_formalism.git diff --git a/kinetic_formalism.Rnw b/kinetic_formalism.Rnw index 8b387e2..fccc4af 100644 --- a/kinetic_formalism.Rnw +++ b/kinetic_formalism.Rnw @@ -59,9 +59,22 @@ \begin{document} %\maketitle -<>= -require(lattice) -require(grid) +<>= +require("lattice") +require("grid") +require("Hmisc") +require("gridBase") +opts_chunk$set(dev="cairo_pdf", + out.width="0.8\\textwidth", + out.height="0.8\\textheight", + out.extra="keepaspectratio") +opts_chunk$set(cache=TRUE, autodep=TRUE) +options(device = function(file, width = 6, height = 6, ...) { + cairo_pdf(tempfile(), width = width, height = height, ...) + }) +to.latex <- function(x){ + gsub("\\\\","\\\\\\\\",latexSN(x)) +} # R in cal / mol K to.kcal <- function(k,temp=300) { gasconst <- 1.985 @@ -72,18 +85,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 $C_{i_\mathrm{ves}}$, -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} - @@ -111,6 +127,7 @@ $\mathrm{n}$, Thus, we have 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$), @@ -148,14 +165,18 @@ $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}}$. -\setkeys{Gin}{width=3.2in} -<>= +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. + + +<>= 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", @@ -191,7 +212,7 @@ 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) @@ -205,7 +226,7 @@ print(wireframe(z~x*y,grid,cuts=50, 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) @@ -262,7 +283,7 @@ of $\Sexpr{format(digits=3,to.kcal(10^(0.13*0.213)))} 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)), @@ -277,7 +298,7 @@ print(wireframe(z~x*y,grid,cuts=50, 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)), @@ -316,14 +337,45 @@ a range of $\Delta \Delta G^\ddagger$ of $\Sexpr{format(digits=3,to.kcal(2^(3.4)))} \frac{\mathrm{kcal}}{\mathrm{mol}}$. - -<>= +While it could be argued that increased length of the monomer could +affect the rate of insertion into the membrane, it's not clear whether +it would increase (by decreasing the number of available hydrogen +bonds, for example) or decrease (by increasing the time taken to fully +insert the acyl chain, for example) the rate of insertion or to what +degree, so we do not take it into account in this formalism. + +\fixme{Incorporate McLean84 here} +From McLean84LIB: Although it is difficult to measure cmc values for +the sparingly soluble lipids used in this study, estimates for +lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds, +1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution; +Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6 +X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X +10-8 M for DMPC was estimated from the linear relationship between ln +cmc and the number of carbons in the PC acyl chain by using data for n += 7, 8, 10, and 16 [summarized in Tanford (1980)]. + +From Nichols85: The magnitude of the dissociation rate constant +decreases by a factor of approximately 3.2 per carbon increase in acyl +chain length (see Table II here) {take into account for the formula; + rz 8/17/2010}. + +From Nichols85: The magnitude of the partition coefficient increases +with acyl chain length [Keq(M-C6-NBD-PC) = (9.8±2.1) X 106 M-1 and Keq +(P-C6-NBD-PC) = (9.4±1.3) x 107 M-1. {take into account for the + formula; rz 8/17/2010}. + +From Nichols85: The association rate constant is independent of acyl +chain length. {take into account for the formula; rz 8/17/2010}. + + +<>= 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", @@ -371,62 +423,6 @@ 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 = 10^{\left(2^{- \left< un_\mathrm{ves} \right> } - -2^{-un_\mathrm{monomer}}\right)^2} - \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(10^((2^-1.7-2^-0)^2)))} -\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with 0 unsaturation -to -$\Sexpr{format(digits=3,to.kcal(10^((2^-1.7-2^-4)^2)))}\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 <- 10^((2^-grid$x-2^-grid$y)^2) -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(10^((2^-grid$x-2^-grid$y)^2)) -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) -@ - -\subsubsection{Unsaturation Backward II} - -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} @@ -441,7 +437,7 @@ $\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-4)^2+1))))}\frac{\mathrm{kc 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))) @@ -453,7 +449,7 @@ print(wireframe(z~x*y,grid,cuts=50, 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)))) @@ -488,7 +484,7 @@ $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) @@ -501,7 +497,7 @@ print(wireframe(z~x*y,grid,cuts=50, 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) @@ -523,9 +519,7 @@ 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. \textcolor{red}{It is this property which the - unsaturation backwards equation does \emph{not} satisfy, which I - think it should.} An equation which satisfies this critera has the +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, @@ -547,7 +541,7 @@ $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(1.3))^2+1))))}\frac{\math 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)) @@ -559,7 +553,7 @@ print(wireframe(z~x*y,grid,cuts=50, 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))) @@ -598,7 +592,7 @@ $\Sexpr{format(digits=3,to.kcal(3.2^abs(24-17.75)))}\frac{\mathrm{kcal}}{\mathrm 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)) @@ -610,7 +604,7 @@ print(wireframe(z~x*y,grid,cuts=50, 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))) @@ -627,6 +621,20 @@ 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} @@ -634,16 +642,18 @@ rm(grid) \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 +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 +\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 length $2$ to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex formation $0$. +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)) @@ -655,7 +665,7 @@ print(wireframe(z~x*y,grid,cuts=50, 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))) @@ -670,6 +680,333 @@ rm(grid) +\subsection{Per-Lipid Kinetic Parameters} + +Each of the 5 lipid types have different kinetic parameters; to the +greatest extent possible, we have derived these from literature. + +\begin{table} + \centering + \begin{tabular}{c c c c c c c} + Type & $k_f$ & $k_b$ & Area (\r{A}$^2$) & Charge & CF1 & Curvature \\ + \hline + PC & $3.7\cdot 10^6$ & $2\cdot 10^{-5}$ & 63 & 0 & 2 & 0.8 \\ + PS & $3.7\cdot 10^6$ & $1.5\cdot 10^{-5}$ & 54 & -1 & 0 & 1 \\ + CHOL & $5.1\cdot 10^7$ & $2.8\cdot 10^{-4}$ & 38 & 0 & -1 & 1.21 \\ + SM & $3.7\cdot 10^6$ & $3.1\cdot 10^{-3}$ & 51 & 0 & 3 & 0.8 \\ + PE & $2.3\cdot 10^6$ & $10^{-5}$ & 55 & 0 & 0 & 1.33 \\ + \end{tabular} + \caption{Kinetic parameters of lipid types} + \label{tab:kinetic_parameters_lipid_types} +\end{table} + +\subsubsection{$k_f$ for lipid types} +For PC, $k_f$ was measured by Nichols85 to be $3.7\cdot 10^6 +\frac{1}{\mathrm{M}\cdot \mathrm{s}}$ by the partitioning of +P-C$_6$-NBD-PC between DOPC vesicles and water. The method utilized by +Nichols85 has the weakness of using NBD-PC, with associated label +perturbations. As similar measures do not exist for SM or PS, we +assume that they have the same $k_f$. For CHOL, Estronca07 found a +value for $k_f$ of $5.1\cdot 10^7 \frac{1}{\mathrm{M}\cdot + \mathrm{s}}$. For PE, Abreu04 found a value for $k_f$ of $2.3\cdot +10^6$. \fixme{I'm missing the notes on these last two papers, so this + isn't correct yet.} + +\subsubsection{$k_b$ for lipid types} + +$k_b$ for PC was measured by Wimley90 using a radioactive label and +large unilammelar vesicles at 30\textdegree C. The other values were +calculated from the experiments of Nichols82 where the ratio of $k_b$ +of different types was measured to that of PC. +See~\fref{tab:kinetic_parameters_lipid_types}. + +assigned accordingly. kb(PS) was assumed to be the same as kb(PG) +given by Nichols82 (also ratio from kb(PC)). kb(SM) is taken from +kb(PC) of Wimley90 (radioactive), and then a ratio of kb(PC)/kb(SM) +taken from Bai97: = 34/2.2 = 15.45; 2.0 x 10-4 x 15.45 = 3.1 x 10-3 s +-1. kb(CHOL) taken from Jones90 (radioactive; POPC LUV; 37°). + +PC 0.89 +PE 0.45 <- from Nichols82 +PG=PS + + +kb PC is from table 2 of Wimley90, where we have a half life of 9.6 +hours for DMPC. \Sexpr{log(2)/(9.6*60*60)}. + + + +\subsubsection{Area for lipid types} + + +From Sampaio05: Besides this work and our own earlier report on the +association of NBD-DMPE with lipid bilayers (Abreu et al., 2004), we +are aware of only one other report in the literature (Nichols, 1985) +in which all the kinetic constants of lipid-derived amphiphile +association with lipid bilayer membranes were experimentally measured. +{indeed; everything is k- !!!; rz} + +From McLean84LIB: Although it is difficult to measure cmc values for +the sparingly soluble lipids used in this study, estimates for +lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds, +1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution; +Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6 +X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X +10-8 M for DMPC was estimated from the linear relationship between ln +cmc and the number of carbons in the PC acyl chain by using data for n += 7, 8, 10, and 16 [summarized in Tanford (1980)]. + +From Toyota08: Recently, several research groups have reported +self-reproducing systems of giant vesicles that undergo a series of +sequential transformations: autonomous growth, self-division, and +chemical reactions to produce membrane constituents within the giant +vesicles.44-47 + +Vesicle sizes of 25 nm for SUV and 150 nm for LUV were mentioned by +Thomas02. + +From Lund-Katz88: Charged and neutral small unilamellar vesicles +composed of either saturated PC, unsaturated PC, or SM had similar +size distributions with diameters of 23 \& 2 nm. + +From Sampaio05LIB: The exchange of lipids and lipid derivatives +between lipid bilayer vesicles has been studied for at least the last +30 years. Most of this work has examined the exchange of amphiphilic +molecules between a donor and an acceptor population. The measured +efflux rates were shown in almost all cases, not surprisingly, to be +first order processes. In all of this work, the insertion rate has +been assumed to be much faster than the efflux rate. Having measured +both the insertion and desorption rate constants for amphiphile +association with membranes, our results show that this assumption is +valid. In several cases reported in the literature, the insertion rate +constant was assumed, although never demonstrated, to be a +diffusion-controlled process. + +(for methods? From McLean84LIB: The activation free energies and free +energies of transfer from self-micelles to water increase by 2.2 and +2.1 kJ mol-' per methylene group, respectively. {see if we can use it + to justify arranging our changed in activating energy around 1 + kcal/mol; rz}). + +Jones90 give diameter of LUV as 100 nm, and of SUV as 20 nm; that +would give the number of molecules per outer leaflet of a vesicle as +1500. + +Form Simard08: Parallel studies with SUV and LUV revealed that +although membrane curvature does have a small effect on the absolute +rates of FA transfer between vesicles, the ΔG of membrane desorption +is unchanged, suggesting that the physical chemical properties which +govern FA desorption are dependent on the dissociating molecule rather +than on membrane curvature. However, disagreements on this fundamental +issue continue (57, 61, 63, 64) + +(methods regarding the curvature effect: Kleinfeld93 showed that the +transfer parameters of long-chain FFA between the lipid vesicles +depend on vesicle curvature and composition. Transfer of stearic acid +is much slower from LUV as compared to SUV). + +From McLean84: In a well-defined experimental system consisting of +unilamellar lipid vesicles, in the absence of protein, the +rate-limiting step for the overall exchange process is desorption +(McLean \& Phillips, 1981). {thus I can take exchange data for the + estimation of k- rz; 8/11/08}. + +\subsubsection{Complex Formation 1} + +From Thomas88a: SM decreases the rate of cholesterol transfer, while +phosphatidylethanolamine (PE) and phosphatidylserine (PS) have no +effect at physiologically significant levels. + + +\section{Simulation Methodology} + +\subsection{Overall Architecture} + +The simulation is currently run by single program written in perl +using various different modules to handle the subsidiary parts. It +produces output for each generation, including the step immediately +preceeding and immediately following a vesicle split (and optionally, +each step) that is written to a state file which contains the state of +the vesicle, environment, kinetic parameters, program invocation +options, time, and various other parameters necessary to recreate the +state vector at a given time. This output file is then read by a +separate program in perl to produce different output which is +generated out-of-band; output which includes graphs and statistical +analysis is performed using R (and various grid graphics modules) +called from the perl program. + +The separation of simulation and output generation allows refining +output, and simulations performed on different versions of the +underlying code to be compared using the same analysis methods and +code. It also allows later simulations to be restarted from a specific +generation, utilizing the same environment. + +Random variables of different distributions are calculated using +Math::Random, which is seeded for each run using entropy from the +linux kernel's urandom device. + +\subsection{Environment Creation} + +\subsubsection{Components} +The environment contains different concentrations of different +components. In the current set of simulations, there are +\Sexpr{1+4*5*7} different components, consisting of PC, PE, PS, SM, +and CHOL, with all lipids except for CHOL having 5 possible +unsaturations rangiong from 0 to 4, and 7 lengths from $12,14,...,22$ +($7\cdot 5\cdot4+1=\Sexpr{1+4*5*7}$). In cases where the environment +has less than the maximum number of components, the components are +selected in order without replacement from a randomly shuffled deck of +component (with the exception of CHOL) represented once until the +desired number of unique components are obtained. CHOL is over +representated $\Sexpr{5*7}$ times to be at the level of other lipid +types, ensures that the probability of CHOL being asbent in the +environment is the same as the probability of one of the other lipid +types (PS, PE, etc.) being entirely absent. This reduces the number of +simulations with a small number of components which are entirely +devoid of CHOL. + +\subsubsection{Concentration} +Once the components of the environment have been selected, the +concentration of thoes components is determined. In experiments where +the environmental concentration is the same across all lipid +components, the concentration is set to $10^{-10}\mathrm{M}$. In other +cases, the environmental concentration is set to a random number from +a gamma distribution with shape parameter 2 and an average of +$10^{-10}\mathrm{M}$. The base concentration ($10^{-10}\mathrm{M}$) +can also be altered in the initialization of the experiment to +specific values for specific lipid types or components. + +\subsection{Initial Vesicle Creation} + +Initial vesicles are seeded by selecting lipid molecules from the +environment until the vesicle reaches a specific starting size. The +vesicle starting size has gamma distribution with shape parameter 2 +and a mean of the experimentally specified starting size, with a +minimum of 5 lipid molecules. Lipid molecules are then selected to be +added to the lipid membrane according to four different methods. In +the constant method, molecules are added in direct proportion to their +concentration in the environment. The uniform method adds molecules in +proportion to their concentration in the environment scaled by a +uniform random value, whereas the random method adds molecules in +proportion to a uniform random value. The final method is a binomial +method, which adjust the porbability of adding a molecule of a +specific component by the concentration of that component, and then +adds the molecules one by one to the membrane. This final method is +also used in the first three methods to add any missing molecules to +the starting vesicle which were unallocated due to the requirement to +add integer numbers of molecules. + +\subsection{Simulation Step} + +Once the environment has been created and the initial vesicle has been +formed, molecules join and leave the vesicle based on the kinetic +parameters and state equation discussed until the vesicle splits +forming two daughter vesicles, one of which the program continues to +follow. + +\subsubsection{Calculation of Vesicle Properties} + +\label{sec:ves_prop} +$S_\mathrm{ves}$ is the surface area of the vesicle, and is the sum of +the surface area of all of the individual lipid components; each lipid +type has a different surface area; we na\"ively assume that the lipid +packing is optimal, and there is no empty space. + +\subsubsection{Joining and Leaving of Lipid Molecules} + +Determining the number of molecules to add to the lipid membrane +($n_i$) requires knowing $k_{fi\mathrm{adj}}$, the surface area of the +vesicle $S_\mathrm{ves}$ (see~\fref{sec:ves_prop}), the time interval +$dt$ during which lipids are added, the base $k_{fi}$, and the +concentration of the monomer in the environment +$\left[C_{i\mathrm{ves}}\right]$ (see~\fref{eq:state}). +$k_{fi\mathrm{adj}}$ is calculated (see~\fref{eq:kf_adj}) based on the +vesicle properties and their hypothesized effect on the rate (in as +many cases as possible, experimentally based) +(see~\fref{sec:forward_adj} for details). $dt$ can be varied +(see~\fref{sec:step_duration}), but for a given step is constant. This +leads to the following: + +$n_i = k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]S_\mathrm{ves}\mathrm{N_A}dt$ + +In the cases where $n_i > 1$, the integer number of molecules is +added. Fractional $n_i$ or the fractional remainder after the addition +of the integer molecules are the probability of adding a specific +molecule, and are compared to a uniformly distributed random value +between 0 and 1. If the random value is less than or equal to the +fractional part of $n_i$, an additional molecule is added. + +Molecules leaving the vesicle are handled in a similar manner, with + +$n_i = k_{bi}k_{bi\mathrm{adj}}C_{i_\mathrm{ves}}\mathrm{N_A}dt$. + +While programatically, the molecule removal happens after the +addition, the properties that each operates on are the same, so they +can be considered to have been added and removed at the same instant. +[This also avoids cases where a removal would have resulted in +negative molecules for a particular type, which is obviously +disallowed.] + +\subsubsection{Step duration} +\label{sec:step_duration} + +$dt$ is the time taken for each step of joining, leaving, and checking +split. In the current implementation, it starts with a value of +$10^{-6}\mathrm{s}$ but this is modified in between each step if the +number of molecules joining or leaving is too large or small. If more +than half of the starting vesicle size molecules join or leave in a +single step, $dt$ is reduced by half. If less than the starting +vesicle size molecules join or leave in 100 steps, $dt$ is doubled. +This is necessary to curtail run times and to automatically adjust to +different experimental runs. + +(In every run seen so far, the initial $dt$ was too small, and was +increased before the first generation occured; at no time was $dt$ too +large.) + +\subsubsection{Vesicle splitting} + +If a vesicle has grown to a size which is more than double the +starting vesicle size, the vesicle splits. More elaborate mechanisms +for determining whether a vesicle should split are of course possible, +but not currently implemented. Molecules are assigned to the daughter +vesicles at random, with each daughter vesicle having an equal chance +of getting a single molecule. The number of molecules to assign to the +first vesicle has binomial distribution with a probability of an event +in each trial of 0.5 and a number of trials equal to the number of +molecules. + +\subsection{Output} + +The environment, initial vesicle, and the state of the vesicle +immediately before and immediately after splitting are stored to disk +to produce later output. + +\section{Analyzing output} + +Analyzing of output is handled by a separate perl program which shares +many common modules with the simulation program. Current output +includes simulation progress, summary tables, summary statistics, and +various graphs. + +\subsection{PCA plots} + +Vesicles have many different axes which contribute to their variation +between subsequent generations; two major groups of axes are the +components and properties of vesicles. Each component in a vesicle is +an axis on its own; it can be measured either as an absolute number of +molecules in each component, or the fraction of molecules of that +component over the total number of molecules; the second approach is +often more convenient, as it allows vesicles of different number of +molecules to be more directly compared (though it hides any affect of +vesicle size). + +In order to visualize the transition of subsequent generations of +vesicles from their initial state in the simulation, to their final +state at the termination of + +\subsection{Carpet plots} + % \bibliographystyle{plainnat}