1 \documentclass[english,12pt]{article}
3 %\usepackage[pdftex]{graphicx}
5 \usepackage[bf]{caption2}
13 \usepackage[light,all]{draftcopy}
15 \usepackage[hyperfigures,backref,bookmarks,colorlinks]{hyperref}
16 \usepackage[sectionbib,sort&compress,square,numbers]{natbib}
17 \usepackage[margin,inline,draft]{fixme}
18 \usepackage[x11names,svgnames]{xcolor}
20 \newenvironment{narrow}[2]{%
22 \setlength{\topsep}{0pt}%
23 \setlength{\leftmargin}{#1}%
24 \setlength{\rightmargin}{#2}%
25 \setlength{\listparindent}{\parindent}%
26 \setlength{\itemindent}{\parindent}%
27 \setlength{\parsep}{\parskip}}%
29 \newenvironment{paperquote}{%
34 \renewcommand{\textfraction}{0.15}
35 \renewcommand{\topfraction}{0.85}
36 \renewcommand{\bottomfraction}{0.65}
37 \renewcommand{\floatpagefraction}{0.60}
38 %\renewcommand{\baselinestretch}{1.8}
39 \newenvironment{enumerate*}%
41 \setlength{\itemsep}{0pt}%
42 \setlength{\parskip}{0pt}}%
44 \newenvironment{itemize*}%
46 \setlength{\itemsep}{0pt}%
47 \setlength{\parskip}{0pt}}%
55 \author{Don Armstrong}
56 \title{OOL Kinetic Formalisms}
62 <<results=hide,echo=FALSE>>=
67 to.latex <- function(x){
68 gsub("\\\\","\\\\\\\\",latexSN(x))
71 to.kcal <- function(k,temp=300) {
73 return(-gasconst*temp*log(k)/1000)
77 \section{State Equation}
78 % double check this with the bits in the paper
80 The base forward kinetic parameter for the $i$th component is
81 $k_{\mathrm{f}i}$ and is dependent on the particular lipid type (PC,
82 PS, SM, etc.). The forward adjustment parameter,
83 $k_{\mathrm{f}i\mathrm{adj}}$, is based on the properties of the
84 vesicle and the specific component (type, length, unsaturation, etc.)
85 (see Equation~\ref{eq:kf_adj}, and
86 Section~\ref{sec:kinetic_adjustments}).
87 $\left[C_{i_\mathrm{monomer}}\right]$ is the molar concentration of
88 monomer of the $i$th component. $\left[S_\mathrm{vesicle}\right]$ is
89 the surface area of the vesicle per volume. The base backwards kinetic
90 parameter for the $i$th component is $k_{\mathrm{b}i}$ and its
91 adjustment parameter $k_{\mathrm{b}i\mathrm{adj}}$ (see
92 Equation~\ref{eq:kb_adj}, and Section~\ref{sec:kinetic_adjustments}).
93 $\left[C_{i_\mathrm{vesicle}}\right]$ is the molar concentration of
94 the $i$th component in the vesicle.
97 \frac{d C_{i_\mathrm{ves}}}{dt} = k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]S_\mathrm{ves} -
98 k_{bi}k_{bi\mathrm{adj}}C_{i_\mathrm{ves}}
102 For $k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]$,
103 $k_{fi}$ has units of $\frac{\mathrm{m}}{\mathrm{s}}$,
104 $k_{fi\mathrm{adj}}$ and $k_{bi\mathrm{adj}}$ are unitless,
105 concentration is in units of $\frac{\mathrm{n}}{\mathrm{L}}$, surface
106 area is in units of $\mathrm{m}^2$, $k_{bi}$ has units of
107 $\frac{1}{\mathrm{s}}$ and $C_{i_\mathrm{ves}}$ has units of
108 $\mathrm{n}$, Thus, we have
111 \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} -
112 \frac{1}{\mathrm{s}} \mathrm{n}
114 \frac{\mathrm{m^3}}{\mathrm{s}} \frac{\mathrm{n}}{\mathrm{L}} \frac{1000\mathrm{L}}{\mathrm{m}^3} - \frac{\mathrm{n}}{\mathrm{s}}=
115 \frac{\mathrm{n}}{\mathrm{s}} = 1000 \frac{\mathrm{n}}{\mathrm{s}} - \frac{\mathrm{n}}{\mathrm{s}}
116 \label{eq:state_units}
119 The 1000 isn't in \fref{eq:state} above, because it is unit-dependent.
121 \subsection{Forward adjustments ($k_{fi\mathrm{adj}}$)}
122 \label{sec:forward_adj}
124 The forward rate constant adjustment, $k_{fi\mathrm{adj}}$ takes into
125 account unsaturation ($un_f$), charge ($ch_f$), curvature ($cu_f$),
126 length ($l_f$), and complex formation ($CF1_f$), each of which are
127 modified depending on the specific specie and the vesicle into which
128 the specie is entering.
131 k_{fi\mathrm{adj}} = un_f \cdot ch_f \cdot cu_f \cdot l_f \cdot CF1_f
136 \subsubsection{Unsaturation Forward}
138 In order for a lipid to be inserted into a membrane, a void has to be
139 formed for it to fill. Voids can be generated by the combination of
140 unsaturated and saturated lipids forming herterogeneous domains. Void
141 formation is increased when the unsaturation of lipids in the vesicle
142 is widely distributed; in other words, the insertion of lipids into
143 the membrane is greater when the standard deviation of the
144 unsaturation is larger. Assuming that an increase in width of the
145 distribution linearly decreases the free energy of activation, the
146 $un_f$ parameter must follow
147 $a^{\mathrm{stdev}\left(un_\mathrm{ves}\right)}$ where $a > 1$, so a
148 convenient starting base for $a$ is $2$:
151 un_f = 2^{\mathrm{stdev}\left(un_\mathrm{ves}\right)}
152 \label{eq:unsaturation_forward}
155 The most common $\mathrm{stdev}\left(un_\mathrm{ves}\right)$ is around
156 $1.5$, which leads to a $\Delta \Delta G^\ddagger$ of
157 $\Sexpr{format(digits=3,to.kcal(2^1.5))}
158 \frac{\mathrm{kcal}}{\mathrm{mol}}$.
160 It is not clear that the unsaturation of the inserted monomer will
161 affect the rate of the insertion positively or negatively, so we do
162 not include a term for it in this formalism.
165 \setkeys{Gin}{width=3.2in}
166 <<fig=TRUE,echo=FALSE,results=hide,width=5,height=5>>=
167 curve(2^x,from=0,to=sd(c(0,4)),
168 main="Unsaturation Forward",
169 xlab="Standard Deviation of Unsaturation of Vesicle",
170 ylab="Unsaturation Forward Adjustment")
172 <<fig=TRUE,echo=FALSE,results=hide,width=5,height=5>>=
173 curve(to.kcal(2^x),from=0,to=sd(c(0,4)),
174 main="Unsaturation forward",
175 xlab="Standard Deviation of Unsaturation of Vesicle",
176 ylab="Unsaturation Forward (kcal/mol)")
181 \subsubsection{Charge Forward}
183 A charged lipid such as PS approaching a vesicle with an average
184 charge of the same sign will experience repulsion, whereas those with
185 different signs will experience attraction, the degree of which is
186 dependent upon the charge of the monomer and the average charge of the
187 vesicle. If either the vesicle or the monomer has no charge, there
188 should be no effect of charge upon the rate. This leads us to the
189 following equation, $a^{-\left<ch_v\right> ch_m}$, where
190 $\left<ch_v\right>$ is the average charge of the vesicle, and $ch_m$
191 is the charge of the monomer. If either $\left<ch_v\right>$ or $ch_m$
192 is 0, the adjustment parameter is 1 (no change), whereas it decreases
193 if both are positive or negative, as the product of two real numbers
194 with the same sign is always positive. A convenient base for $a$ is
195 60, resulting in the following equation:
199 ch_f = 60^{-\left<{ch}_v\right> {ch}_m}
200 \label{eq:charge_forward}
203 The most common $\left<{ch}_v\right>$ is around $-0.165$, which leads to
204 a range of $\Delta \Delta G^\ddagger$ from
205 $\Sexpr{format(digits=3,to.kcal(60^(-.165*-1)))}
206 \frac{\mathrm{kcal}}{\mathrm{mol}}$ to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$.
208 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
209 x <- seq(-1,0,length.out=20)
210 y <- seq(-1,0,length.out=20)
211 grid <- expand.grid(x=x,y=y)
212 grid$z <- as.vector(60^(-outer(x,y)))
213 print(wireframe(z~x*y,grid,cuts=50,
215 scales=list(arrows=FALSE),
216 main="Charge Forward",
217 xlab=list("Average Vesicle Charge",rot=30),
218 ylab=list("Component Charge",rot=-35),
219 zlab=list("Charge Forward",rot=93)))
222 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
223 x <- seq(-1,0,length.out=20)
224 y <- seq(-1,0,length.out=20)
225 grid <- expand.grid(x=x,y=y)
226 grid$z <- as.vector(to.kcal(60^(-outer(x,y))))
227 print(wireframe(z~x*y,grid,cuts=50,
229 scales=list(arrows=FALSE),
230 main="Charge Forward (kcal/mol)",
231 xlab=list("Average Vesicle Charge",rot=30),
232 ylab=list("Component Charge",rot=-35),
233 zlab=list("Charge Forward (kcal/mol)",rot=93)))
239 \subsubsection{Curvature Forward}
241 Curvature is a measure of the intrinsic propensity of specific lipids
242 to form micelles (positive curvature), inverted micelles (negative
243 curvature), or planar sheets (zero curvature). In this formalism,
244 curvature is measured as the ratio of the size of the head to that of
245 the base, so negative curvature is bounded by $(0,1)$, zero curvature
246 is 1, and positive curvature is bounded by $(1,\infty)$. The curvature
247 can be transformed into the typical postive/negative mapping using
248 $\log$, which has the additional property of making the range of
249 positive and negative curvature equal, and distributed about 0.
251 As in the case of unsaturation, void formation is increased by the
252 presence of lipids with mismatched curvature. Thus, a larger
253 distribution of curvature in the vesicle increases the rate of lipid
254 insertion into the vesicle. However, a species with curvature $e^{-1}$
255 will cancel out a species with curvature $e$, so we have to log
256 transform (turning these into -1 and 1), then take the absolute value
257 (1 and 1), and finally measure the width of the distribution. Thus, by
258 using the log transform to make the range of the lipid curvature equal
259 between positive and negative, and taking the average to cancel out
260 exactly mismatched curvatures, we come to an equation with the shape
261 $a^{\left<\log cu_\mathrm{vesicle}\right>}$. A convenient base for $a$
266 % cu_f = 10^{\mathrm{stdev}\left|\log cu_\mathrm{vesicle}\right|}
267 cu_f = 10^{\left|\left<\log cu_\mathrm{vesicle} \right>\right|\mathrm{stdev} \left|\log cu_\mathrm{vesicle}\right|}
268 \label{eq:curvature_forward}
271 The most common $\left|\left<\log {cu}_v\right>\right|$ is around
272 $0.013$, which with the most common $\mathrm{stdev} \log
273 cu_\mathrm{vesicle}$ of $0.213$ leads to a $\Delta \Delta G^\ddagger$
274 of $\Sexpr{format(digits=3,to.kcal(10^(0.13*0.213)))}
275 \frac{\mathrm{kcal}}{\mathrm{mol}}$. This is a consequence of the
276 relatively matched curvatures in our environment.
278 % 1.5 to 0.75 3 to 0.33
279 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
280 grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
281 y sd(abs(log(c(1,0.33)))),sd(abs(log(c(0.33,3)))))),length.out=20),
282 y=seq(0,max(c(mean(log(c(1,3)),
283 mean(log(c(1,0.33))),
284 mean(log(c(0.33,3)))))),length.out=20))
285 grid$z <- 10^(grid$x*grid$y)
286 print(wireframe(z~x*y,grid,cuts=50,
288 scales=list(arrows=FALSE),
289 xlab=list("Vesicle stdev log curvature",rot=30),
290 ylab=list("Vesicle average log curvature",rot=-35),
291 zlab=list("Vesicle Curvature Forward",rot=93)))
294 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
295 grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
296 sd(abs(log(c(1,0.33)))),sd(abs(log(c(0.33,3)))))),length.out=20),
297 y=seq(0,max(c(mean(log(c(1,3)),
298 mean(log(c(1,0.33))),
299 mean(log(c(0.33,3)))))),length.out=20))
300 grid$z <- to.kcal(10^(grid$x*grid$y))
301 print(wireframe(z~x*y,grid,cuts=50,
303 scales=list(arrows=FALSE),
304 xlab=list("Vesicle stdev log curvature",rot=30),
305 ylab=list("Vesicle average log curvature",rot=-35),
306 zlab=list("Vesicle Curvature Forward (kcal/mol)",rot=93)))
311 \subsubsection{Length Forward}
313 As in the case of unsaturation, void formation is easier when vesicles
314 are made up of components of widely different lengths. Thus, when the
315 width of the distribution of lengths is larger, the forward rate
316 should be greater as well, leading us to an equation of the form
317 $x^{\mathrm{stdev} l_\mathrm{ves}}$, where $\mathrm{stdev}
318 l_\mathrm{ves}$ is the standard deviation of the length of the
319 components of the vesicle, which has a maximum possible value of 8 and
320 a minimum of 0 in this set of experiments. A convenient base for $x$
324 l_f = 2^{\mathrm{stdev} l_\mathrm{ves}}
325 \label{eq:length_forward}
328 The most common $\mathrm{stdev} l_\mathrm{ves}$ is around $3.4$, which leads to
329 a range of $\Delta \Delta G^\ddagger$ of
330 $\Sexpr{format(digits=3,to.kcal(2^(3.4)))}
331 \frac{\mathrm{kcal}}{\mathrm{mol}}$.
333 While it could be argued that increased length of the monomer could
334 affect the rate of insertion into the membrane, it's not clear whether
335 it would increase (by decreasing the number of available hydrogen
336 bonds, for example) or decrease (by increasing the time taken to fully
337 insert the acyl chain, for example) the rate of insertion or to what
338 degree, so we do not take it into account in this formalism.
340 \fixme{Incorporate McLean84 here}
341 From McLean84LIB: Although it is difficult to measure cmc values for
342 the sparingly soluble lipids used in this study, estimates for
343 lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds,
344 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution;
345 Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6
346 X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X
347 10-8 M for DMPC was estimated from the linear relationship between ln
348 cmc and the number of carbons in the PC acyl chain by using data for n
349 = 7, 8, 10, and 16 [summarized in Tanford (1980)].
351 From Nichols85: The magnitude of the dissociation rate constant
352 decreases by a factor of approximately 3.2 per carbon increase in acyl
353 chain length (see Table II here) {take into account for the formula;
356 From Nichols85: The magnitude of the partition coefficient increases
357 with acyl chain length [Keq(M-C6-NBD-PC) = (9.8±2.1) X 106 M-1 and Keq
358 (P-C6-NBD-PC) = (9.4±1.3) x 107 M-1. {take into account for the
359 formula; rz 8/17/2010}.
361 From Nichols85: The association rate constant is independent of acyl
362 chain length. {take into account for the formula; rz 8/17/2010}.
365 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=5>>=
366 curve(2^x,from=0,to=sd(c(12,24)),
367 main="Length forward",
368 xlab="Standard Deviation of Length of Vesicle",
369 ylab="Length Forward Adjustment")
371 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=5>>=
372 curve(to.kcal(2^x),from=0,to=sd(c(12,24)),
373 main="Length forward",
374 xlab="Standard Deviation of Length of Vesicle",
375 ylab="Length Forward Adjustment (kcal/mol)")
379 \subsubsection{Complex Formation}
380 There is no contribution of complex formation to the forward reaction
381 rate in the current formalism.
385 \label{eq:complex_formation_forward}
388 \subsection{Backward adjustments ($k_{bi\mathrm{adj}}$)}
390 Just as the forward rate constant adjustment $k_{fi\mathrm{adj}}$
391 does, the backwards rate constant adjustment $k_{bi\mathrm{adj}}$
392 takes into account unsaturation ($un_b$), charge ($ch_b$), curvature
393 ($cu_b$), length ($l_b$), and complex formation ($CF1_b$), each of
394 which are modified depending on the specific specie and the vesicle
395 into which the specie is entering:
399 k_{bi\mathrm{adj}} = un_b \cdot ch_b \cdot cu_b \cdot l_b \cdot CF1_b
403 \subsubsection{Unsaturation Backward}
405 Unsaturation also influences the ability of a lipid molecule to leave
406 a membrane. If a molecule has an unsaturation level which is different
407 from the surrounding membrane, it will be more likely to leave the
408 membrane. The more different the unsaturation level is, the greater
409 the propensity for the lipid molecule to leave. However, a vesicle
410 with some unsaturation is more favorable for lipids with more
411 unsaturation than the equivalent amount of less unsatuturation, so the
412 difference in energy between unsaturation is not linear. Therefore, an
413 equation with the shape
414 $x^{\left| y^{-\left< un_\mathrm{ves}\right> }-y^{-un_\mathrm{monomer}} \right| }$
415 where $\left<un_\mathrm{ves}\right>$ is the average unsaturation of
416 the vesicle, and $un_\mathrm{monomer}$ is the average unsaturation. In
417 this equation, as the average unsaturation of the vesicle is larger,
420 un_b = 7^{1-\left(20\left(2^{-\left<un_\mathrm{vesicle} \right>} - 2^{-un_\mathrm{monomer}} \right)^2+1\right)^{-1}}
421 \label{eq:unsaturation_backward}
424 The most common $\left<un_\mathrm{ves}\right>$ is around $1.7$, which leads to
425 a range of $\Delta \Delta G^\ddagger$ from
426 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-0)^2+1))))}
427 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with 0 unsaturation
429 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-4)^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
430 for monomers with 4 unsaturations.
433 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
434 grid <- expand.grid(x=seq(0,4,length.out=20),
435 y=seq(0,4,length.out=20))
436 grid$z <- (7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1)))
437 print(wireframe(z~x*y,grid,cuts=50,
439 scales=list(arrows=FALSE),
440 xlab=list("Average Vesicle Unsaturation",rot=30),
441 ylab=list("Monomer Unsaturation",rot=-35),
442 zlab=list("Unsaturation Backward",rot=93)))
445 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
446 grid <- expand.grid(x=seq(0,4,length.out=20),
447 y=seq(0,4,length.out=20))
448 grid$z <- to.kcal((7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1))))
449 print(wireframe(z~x*y,grid,cuts=50,
451 scales=list(arrows=FALSE),
452 xlab=list("Average Vesicle Unsaturation",rot=30),
453 ylab=list("Monomer Unsaturation",rot=-35),
454 zlab=list("Unsaturation Backward (kcal/mol)",rot=93)))
461 \subsubsection{Charge Backwards}
462 As in the case of monomers entering a vesicle, monomers leaving a
463 vesicle leave faster if their charge has the same sign as the average
464 charge vesicle. An equation of the form $ch_b = a^{\left<ch_v\right>
465 ch_m}$ is then appropriate, and using a base of $a=20$ yields:
468 ch_b = 20^{\left<{ch}_v\right> {ch}_m}
469 \label{eq:charge_backwards}
472 The most common $\left<ch_v\right>$ is around $-0.164$, which leads to
473 a range of $\Delta \Delta G^\ddagger$ from
474 $\Sexpr{format(digits=3,to.kcal(20^(-.164*-1)))}
475 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with charge $-1$ to
476 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$
477 for monomers with charge $0$.
480 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
481 x <- seq(-1,0,length.out=20)
482 y <- seq(-1,0,length.out=20)
483 grid <- expand.grid(x=x,y=y)
484 grid$z <- as.vector(20^(outer(x,y)))
485 print(wireframe(z~x*y,grid,cuts=50,
487 scales=list(arrows=FALSE),
488 xlab=list("Average Vesicle Charge",rot=30),
489 ylab=list("Component Charge",rot=-35),
490 zlab=list("Charge Backwards",rot=93)))
493 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
494 x <- seq(-1,0,length.out=20)
495 y <- seq(-1,0,length.out=20)
496 grid <- expand.grid(x=x,y=y)
497 grid$z <- to.kcal(as.vector(20^(outer(x,y))))
498 print(wireframe(z~x*y,grid,cuts=50,
500 scales=list(arrows=FALSE),
501 xlab=list("Average Vesicle Charge",rot=30),
502 ylab=list("Component Charge",rot=-35),
503 zlab=list("Charge Backwards (kcal/mol)",rot=93)))
508 \subsubsection{Curvature Backwards}
510 The less a monomer's intrinsic curvature matches the average curvature
511 of the vesicle in which it is in, the greater its rate of efflux. If
512 the difference is 0, $cu_f$ needs to be one. To map negative and
513 positive curvature to the same range, we also need take the logarithm.
514 Increasing mismatches in curvature increase the rate of efflux, but
515 asymptotically. An equation which satisfies this critera has the
516 form $cu_f = a^{1-\left(b\left( \left< \log cu_\mathrm{vesicle} \right>
517 -\log cu_\mathrm{monomer}\right)^2+1\right)^{-1}}$. An
518 alternative form would use the aboslute value of the difference,
519 however, this yields a cusp and sharp increase about the curvature
520 equilibrium, which is decidedly non-elegant. We have chosen bases of
524 cu_b = 7^{1-\left(20\left(\left< \log cu_\mathrm{vesicle} \right> -\log cu_\mathrm{monomer} \right)^2+1\right)^{-1}}
525 \label{eq:curvature_backwards}
528 The most common $\left<\log cu_\mathrm{ves}\right>$ is around $-0.013$, which leads to
529 a range of $\Delta \Delta G^\ddagger$ from
530 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(0.8))^2+1))))}
531 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature 0.8
533 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(1.3))^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
534 for monomers with curvature 1.3 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature near 1.
537 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
538 grid <- expand.grid(x=seq(0.8,1.33,length.out=20),
539 y=seq(0.8,1.33,length.out=20))
540 grid$z <- 7^(1-1/(20*(log(grid$x)-log(grid$y))^2+1))
541 print(wireframe(z~x*y,grid,cuts=50,
543 scales=list(arrows=FALSE),
544 xlab=list("Vesicle Curvature",rot=30),
545 ylab=list("Monomer Curvature",rot=-35),
546 zlab=list("Curvature Backward",rot=93)))
549 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
550 grid <- expand.grid(x=seq(0.8,1.33,length.out=20),
551 y=seq(0.8,1.33,length.out=20))
552 grid$z <- to.kcal(7^(1-1/(20*(log(grid$x)-log(grid$y))^2+1)))
553 print(wireframe(z~x*y,grid,cuts=50,
555 scales=list(arrows=FALSE),
556 xlab=list("Vesicle Curvature",rot=30),
557 ylab=list("Monomer Curvature",rot=-35),
558 zlab=list("Curvature Backward (kcal/mol)",rot=93)))
563 \subsubsection{Length Backwards}
565 In a model membrane, the dissociation constant increases by a factor
566 of approximately 3.2 per carbon decrease in acyl chain length (Nichols
567 1985). Unfortunatly, the experimental data known to us only measures
568 chain length less than or equal to the bulk lipid, and does not exceed
569 it, and is only known for one bulk lipid species (DOPC). We assume
570 then, that the increase is in relationship to the average vesicle, and
571 that lipids with larger acyl chain length will also show an increase
572 in their dissociation constant.
575 l_b = 3.2^{\left|\left<l_\mathrm{ves}\right>-l_\mathrm{monomer}\right|}
576 \label{eq:length_backward}
579 The most common $\left<\log l_\mathrm{ves}\right>$ is around $17.75$, which leads to
580 a range of $\Delta \Delta G^\ddagger$ from
581 $\Sexpr{format(digits=3,to.kcal(3.2^abs(12-17.75)))}
582 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with length 12
584 $\Sexpr{format(digits=3,to.kcal(3.2^abs(24-17.75)))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
585 for monomers with length 24 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature near 18.
588 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
589 grid <- expand.grid(x=seq(12,24,length.out=20),
590 y=seq(12,24,length.out=20))
591 grid$z <- 3.2^(abs(grid$x-grid$y))
592 print(wireframe(z~x*y,grid,cuts=50,
594 scales=list(arrows=FALSE),
595 xlab=list("Average Vesicle Length",rot=30),
596 ylab=list("Monomer Length",rot=-35),
597 zlab=list("Length Backward",rot=93)))
600 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
601 grid <- expand.grid(x=seq(12,24,length.out=20),
602 y=seq(12,24,length.out=20))
603 grid$z <- to.kcal(3.2^(abs(grid$x-grid$y)))
604 print(wireframe(z~x*y,grid,cuts=50,
606 scales=list(arrows=FALSE),
607 xlab=list("Average Vesicle Length",rot=30),
608 ylab=list("Monomer Length",rot=-35),
609 zlab=list("Length Backward (kcal/mol)",rot=93)))
615 \subsubsection{Complex Formation Backward}
617 Complex formation describes the interaction between CHOL and PC or SM,
618 where PC or SM protects the hydroxyl group of CHOL from interactions
619 with water, the ``Umbrella Model''. PC ($CF1=2$) can interact with two
620 CHOL, and SM ($CF1=3$) with three CHOL ($CF1=-1$). If the average of
621 $CF1$ is positive (excess of SM and PC with regards to complex
622 formation), species with negative $CF1$ (CHOL) will be retained. If
623 average $CF1$ is negative, species with positive $CF1$ are retained.
624 An equation which has this property is
625 $CF1_b=a^{\left<CF1_\mathrm{ves}\right>
626 CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{ves}\right>
627 CF1_\mathrm{monomer}\right|}$, where difference of the exponent is
628 zero if the average $CF1$ and the $CF1$ of the specie have the same
629 sign, or double the product if the signs are different. A convenient
630 base for $a$ is $1.5$.
634 CF1_b=1.5^{\left<CF1_\mathrm{ves}\right> CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{ves}\right> CF1_\mathrm{monomer}\right|}
635 \label{eq:complex_formation_backward}
638 The most common $\left<CF1_\mathrm{ves}\right>$ is around $0.925$,
639 which leads to a range of $\Delta \Delta G^\ddagger$ from
640 $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*-1-abs(0.925*-1))))}
641 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
643 $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*2-abs(0.925*2))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
644 for monomers with complex formation $2$ to
645 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
649 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
650 grid <- expand.grid(x=seq(-1,3,length.out=20),
651 y=seq(-1,3,length.out=20))
652 grid$z <- 1.5^(grid$x*grid$y-abs(grid$x*grid$y))
653 print(wireframe(z~x*y,grid,cuts=50,
655 scales=list(arrows=FALSE),
656 xlab=list("Vesicle Complex Formation",rot=30),
657 ylab=list("Monomer Complex Formation",rot=-35),
658 zlab=list("Complex Formation Backward",rot=93)))
661 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
662 grid <- expand.grid(x=seq(-1,3,length.out=20),
663 y=seq(-1,3,length.out=20))
664 grid$z <- to.kcal(1.5^(grid$x*grid$y-abs(grid$x*grid$y)))
665 print(wireframe(z~x*y,grid,cuts=50,
667 scales=list(arrows=FALSE),
668 xlab=list("Vesicle Complex Formation",rot=30),
669 ylab=list("Monomer Complex Formation",rot=-35),
670 zlab=list("Complex Formation Backward (kcal/mol)",rot=93)))
676 \subsection{Per-Lipid Kinetic Parameters}
678 Each of the 5 lipid types have different kinetic parameters; to the
679 greatest extent possible, we have derived these from literature.
683 \begin{tabular}{c c c c c c c}
684 Type & $k_f$ & $k_b$ & Area (\r{A}$^2$) & Charge & CF1 & Curvature \\
686 PC & $3.7\cdot 10^6$ & $2\cdot 10^{-5}$ & 63 & 0 & 2 & 0.8 \\
687 PS & $3.7\cdot 10^6$ & $1.5\cdot 10^{-5}$ & 54 & -1 & 0 & 1 \\
688 CHOL & $5.1\cdot 10^7$ & $2.8\cdot 10^{-4}$ & 38 & 0 & -1 & 1.21 \\
689 SM & $3.7\cdot 10^6$ & $3.1\cdot 10^{-3}$ & 51 & 0 & 3 & 0.8 \\
690 PE & $2.3\cdot 10^6$ & $10^{-5}$ & 55 & 0 & 0 & 1.33 \\
692 \caption{Kinetic parameters of lipid types}
693 \label{tab:kinetic_parameters_lipid_types}
696 \subsubsection{$k_f$ for lipid types}
697 For PC, $k_f$ was measured by Nichols85 to be $3.7\cdot 10^6
698 \frac{1}{\mathrm{M}\cdot \mathrm{s}}$ by the partitioning of
699 P-C$_6$-NBD-PC between DOPC vesicles and water. The method utilized by
700 Nichols85 has the weakness of using NBD-PC, with associated label
701 perturbations. As similar measures do not exist for SM or PS, we
702 assume that they have the same $k_f$. For CHOL, Estronca07 found a
703 value for $k_f$ of $5.1\cdot 10^7 \frac{1}{\mathrm{M}\cdot
704 \mathrm{s}}$. For PE, Abreu04 found a value for $k_f$ of $2.3\cdot
705 10^6$. \fixme{I'm missing the notes on these last two papers, so this
708 \subsubsection{$k_b$ for lipid types}
710 $k_b$ for PC was measured by Wimley90 using a radioactive label and
711 large unilammelar vesicles at 30\textdegree C. The other values were
712 calculated from the experiments of Nichols82 where the ratio of $k_b$
713 of different types was measured to that of PC.
714 See~\fref{tab:kinetic_parameters_lipid_types}.
716 assigned accordingly. kb(PS) was assumed to be the same as kb(PG)
717 given by Nichols82 (also ratio from kb(PC)). kb(SM) is taken from
718 kb(PC) of Wimley90 (radioactive), and then a ratio of kb(PC)/kb(SM)
719 taken from Bai97: = 34/2.2 = 15.45; 2.0 x 10-4 x 15.45 = 3.1 x 10-3 s
720 -1. kb(CHOL) taken from Jones90 (radioactive; POPC LUV; 37°).
723 PE 0.45 <- from Nichols82
727 kb PC is from table 2 of Wimley90, where we have a half life of 9.6
728 hours for DMPC. \Sexpr{log(2)/(9.6*60*60)}.
732 \subsubsection{Area for lipid types}
735 From Sampaio05: Besides this work and our own earlier report on the
736 association of NBD-DMPE with lipid bilayers (Abreu et al., 2004), we
737 are aware of only one other report in the literature (Nichols, 1985)
738 in which all the kinetic constants of lipid-derived amphiphile
739 association with lipid bilayer membranes were experimentally measured.
740 {indeed; everything is k- !!!; rz}
742 From McLean84LIB: Although it is difficult to measure cmc values for
743 the sparingly soluble lipids used in this study, estimates for
744 lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds,
745 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution;
746 Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6
747 X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X
748 10-8 M for DMPC was estimated from the linear relationship between ln
749 cmc and the number of carbons in the PC acyl chain by using data for n
750 = 7, 8, 10, and 16 [summarized in Tanford (1980)].
752 From Toyota08: Recently, several research groups have reported
753 self-reproducing systems of giant vesicles that undergo a series of
754 sequential transformations: autonomous growth, self-division, and
755 chemical reactions to produce membrane constituents within the giant
758 Vesicle sizes of 25 nm for SUV and 150 nm for LUV were mentioned by
761 From Lund-Katz88: Charged and neutral small unilamellar vesicles
762 composed of either saturated PC, unsaturated PC, or SM had similar
763 size distributions with diameters of 23 \& 2 nm.
765 From Sampaio05LIB: The exchange of lipids and lipid derivatives
766 between lipid bilayer vesicles has been studied for at least the last
767 30 years. Most of this work has examined the exchange of amphiphilic
768 molecules between a donor and an acceptor population. The measured
769 efflux rates were shown in almost all cases, not surprisingly, to be
770 first order processes. In all of this work, the insertion rate has
771 been assumed to be much faster than the efflux rate. Having measured
772 both the insertion and desorption rate constants for amphiphile
773 association with membranes, our results show that this assumption is
774 valid. In several cases reported in the literature, the insertion rate
775 constant was assumed, although never demonstrated, to be a
776 diffusion-controlled process.
778 (for methods? From McLean84LIB: The activation free energies and free
779 energies of transfer from self-micelles to water increase by 2.2 and
780 2.1 kJ mol-' per methylene group, respectively. {see if we can use it
781 to justify arranging our changed in activating energy around 1
784 Jones90 give diameter of LUV as 100 nm, and of SUV as 20 nm; that
785 would give the number of molecules per outer leaflet of a vesicle as
788 Form Simard08: Parallel studies with SUV and LUV revealed that
789 although membrane curvature does have a small effect on the absolute
790 rates of FA transfer between vesicles, the ΔG of membrane desorption
791 is unchanged, suggesting that the physical chemical properties which
792 govern FA desorption are dependent on the dissociating molecule rather
793 than on membrane curvature. However, disagreements on this fundamental
794 issue continue (57, 61, 63, 64)
796 (methods regarding the curvature effect: Kleinfeld93 showed that the
797 transfer parameters of long-chain FFA between the lipid vesicles
798 depend on vesicle curvature and composition. Transfer of stearic acid
799 is much slower from LUV as compared to SUV).
801 From McLean84: In a well-defined experimental system consisting of
802 unilamellar lipid vesicles, in the absence of protein, the
803 rate-limiting step for the overall exchange process is desorption
804 (McLean \& Phillips, 1981). {thus I can take exchange data for the
805 estimation of k- rz; 8/11/08}.
807 \subsubsection{Complex Formation 1}
809 From Thomas88a: SM decreases the rate of cholesterol transfer, while
810 phosphatidylethanolamine (PE) and phosphatidylserine (PS) have no
811 effect at physiologically significant levels.
814 \section{Simulation Methodology}
816 \subsection{Overall Architecture}
818 The simulation is currently run by single program written in perl
819 using various different modules to handle the subsidiary parts. It
820 produces output for each generation, including the step immediately
821 preceeding and immediately following a vesicle split (and optionally,
822 each step) that is written to a state file which contains the state of
823 the vesicle, environment, kinetic parameters, program invocation
824 options, time, and various other parameters necessary to recreate the
825 state vector at a given time. This output file is then read by a
826 separate program in perl to produce different output which is
827 generated out-of-band; output which includes graphs and statistical
828 analysis is performed using R (and various grid graphics modules)
829 called from the perl program.
831 The separation of simulation and output generation allows refining
832 output, and simulations performed on different versions of the
833 underlying code to be compared using the same analysis methods and
834 code. It also allows later simulations to be restarted from a specific
835 generation, utilizing the same environment.
837 Random variables of different distributions are calculated using
838 Math::Random, which is seeded for each run using entropy from the
839 linux kernel's urandom device.
841 \subsection{Environment Creation}
843 \subsubsection{Components}
844 The environment contains different concentrations of different
845 components. In the current set of simulations, there are
846 \Sexpr{1+4*5*7} different components, consisting of PC, PE, PS, SM,
847 and CHOL, with all lipids except for CHOL having 5 possible
848 unsaturations rangiong from 0 to 4, and 7 lengths from $12,14,...,22$
849 ($7\cdot 5\cdot4+1=\Sexpr{1+4*5*7}$). In cases where the environment
850 has less than the maximum number of components, the components are
851 selected in order without replacement from a randomly shuffled deck of
852 component (with the exception of CHOL) represented once until the
853 desired number of unique components are obtained. CHOL is over
854 representated $\Sexpr{5*7}$ times to be at the level of other lipid
855 types, ensures that the probability of CHOL being asbent in the
856 environment is the same as the probability of one of the other lipid
857 types (PS, PE, etc.) being entirely absent. This reduces the number of
858 simulations with a small number of components which are entirely
861 \subsubsection{Concentration}
862 Once the components of the environment have been selected, the
863 concentration of thoes components is determined. In experiments where
864 the environmental concentration is the same across all lipid
865 components, the concentration is set to $10^{-10}\mathrm{M}$. In other
866 cases, the environmental concentration is set to a random number from
867 a gamma distribution with shape parameter 2 and an average of
868 $10^{-10}\mathrm{M}$. The base concentration ($10^{-10}\mathrm{M}$)
869 can also be altered in the initialization of the experiment to
870 specific values for specific lipid types or components.
872 \subsection{Initial Vesicle Creation}
874 Initial vesicles are seeded by selecting lipid molecules from the
875 environment until the vesicle reaches a specific starting size. The
876 vesicle starting size has gamma distribution with shape parameter 2
877 and a mean of the experimentally specified starting size, with a
878 minimum of 5 lipid molecules. Lipid molecules are then selected to be
879 added to the lipid membrane according to four different methods. In
880 the constant method, molecules are added in direct proportion to their
881 concentration in the environment. The uniform method adds molecules in
882 proportion to their concentration in the environment scaled by a
883 uniform random value, whereas the random method adds molecules in
884 proportion to a uniform random value. The final method is a binomial
885 method, which adjust the porbability of adding a molecule of a
886 specific component by the concentration of that component, and then
887 adds the molecules one by one to the membrane. This final method is
888 also used in the first three methods to add any missing molecules to
889 the starting vesicle which were unallocated due to the requirement to
890 add integer numbers of molecules.
892 \subsection{Simulation Step}
894 Once the environment has been created and the initial vesicle has been
895 formed, molecules join and leave the vesicle based on the kinetic
896 parameters and state equation discussed until the vesicle splits
897 forming two daughter vesicles, one of which the program continues to
900 \subsubsection{Calculation of Vesicle Properties}
903 $S_\mathrm{ves}$ is the surface area of the vesicle, and is the sum of
904 the surface area of all of the individual lipid components; each lipid
905 type has a different surface area; we na\"ively assume that the lipid
906 packing is optimal, and there is no empty space.
908 \subsubsection{Joining and Leaving of Lipid Molecules}
910 Determining the number of molecules to add to the lipid membrane
911 ($n_i$) requires knowing $k_{fi\mathrm{adj}}$, the surface area of the
912 vesicle $S_\mathrm{ves}$ (see~\fref{sec:ves_prop}), the time interval
913 $dt$ during which lipids are added, the base $k_{fi}$, and the
914 concentration of the monomer in the environment
915 $\left[C_{i\mathrm{ves}}\right]$ (see~\fref{eq:state}).
916 $k_{fi\mathrm{adj}}$ is calculated (see~\fref{eq:kf_adj}) based on the
917 vesicle properties and their hypothesized effect on the rate (in as
918 many cases as possible, experimentally based)
919 (see~\fref{sec:forward_adj} for details). $dt$ can be varied
920 (see~\fref{sec:step_duration}), but for a given step is constant. This
921 leads to the following:
923 $n_i = k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]S_\mathrm{ves}\mathrm{N_A}dt$
925 In the cases where $n_i > 1$, the integer number of molecules is
926 added. Fractional $n_i$ or the fractional remainder after the addition
927 of the integer molecules are the probability of adding a specific
928 molecule, and are compared to a uniformly distributed random value
929 between 0 and 1. If the random value is less than or equal to the
930 fractional part of $n_i$, an additional molecule is added.
932 Molecules leaving the vesicle are handled in a similar manner, with
934 $n_i = k_{bi}k_{bi\mathrm{adj}}C_{i_\mathrm{ves}}\mathrm{N_A}dt$.
936 While programatically, the molecule removal happens after the
937 addition, the properties that each operates on are the same, so they
938 can be considered to have been added and removed at the same instant.
939 [This also avoids cases where a removal would have resulted in
940 negative molecules for a particular type, which is obviously
943 \subsubsection{Step duration}
944 \label{sec:step_duration}
946 $dt$ is the time taken for each step of joining, leaving, and checking
947 split. In the current implementation, it starts with a value of
948 $10^{-6}\mathrm{s}$ but this is modified in between each step if the
949 number of molecules joining or leaving is too large or small. If more
950 than half of the starting vesicle size molecules join or leave in a
951 single step, $dt$ is reduced by half. If less than the starting
952 vesicle size molecules join or leave in 100 steps, $dt$ is doubled.
953 This is necessary to curtail run times and to automatically adjust to
954 different experimental runs.
956 (In every run seen so far, the initial $dt$ was too small, and was
957 increased before the first generation occured; at no time was $dt$ too
960 \subsubsection{Vesicle splitting}
962 If a vesicle has grown to a size which is more than double the
963 starting vesicle size, the vesicle splits. More elaborate mechanisms
964 for determining whether a vesicle should split are of course possible,
965 but not currently implemented. Molecules are assigned to the daughter
966 vesicles at random, with each daughter vesicle having an equal chance
967 of getting a single molecule. The number of molecules to assign to the
968 first vesicle has binomial distribution with a probability of an event
969 in each trial of 0.5 and a number of trials equal to the number of
974 The environment, initial vesicle, and the state of the vesicle
975 immediately before and immediately after splitting are stored to disk
976 to produce later output.
978 \section{Analyzing output}
980 Analyzing of output is handled by a separate perl program which shares
981 many common modules with the simulation program. Current output
982 includes simulation progress, summary tables, summary statistics, and
985 \subsection{PCA plots}
987 Vesicles have many different axes which contribute to their variation
988 between subsequent generations; two major groups of axes are the
989 components and properties of vesicles. Each component in a vesicle is
990 an axis on its own; it can be measured either as an absolute number of
991 molecules in each component, or the fraction of molecules of that
992 component over the total number of molecules; the second approach is
993 often more convenient, as it allows vesicles of different number of
994 molecules to be more directly compared (though it hides any affect of
997 In order to visualize the transition of subsequent generations of
998 vesicles from their initial state in the simulation, to their final
999 state at the termination of
1001 \subsection{Carpet plots}
1005 % \bibliographystyle{plainnat}
1006 % \bibliography{references.bib}