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>>=
66 to.kcal <- function(k,temp=300) {
68 return(-gasconst*temp*log(k)/1000)
72 \section{State Equation}
73 % double check this with the bits in the paper
75 Given a base forward kinetic parameter for the $i$th specie $k_{fi}$
76 (which is dependent on lipid type, that is PC, PE, PS, etc.), an
77 adjustment parameter $k_{fi\mathrm{adj}}$ based on the vesicle and the
78 specific specie (length, unsaturation, etc.) (see~\fref{eq:kf_adj}),
79 the molar concentration of monomer of the $i$th specie
80 $\left[C_{i_\mathrm{monomer}}\right]$, the surface area of the vesicle
81 $S_\mathrm{ves}$, the base backwards kinetic parameter for the $i$th
82 specie $k_{bi}$ which is also dependent on lipid type, its adjustment
83 parameter $k_{bi\mathrm{adj}}$ (see~\fref{eq:kb_adj}), and the molar
84 concentration of the $i$th specie in the vesicle
85 $\left[C_{i_\mathrm{ves}}\right]$, the change in concentration of the
86 $i$th specie in the vesicle per change in time $\frac{d
87 C_{i_\mathrm{ves}}}{dt}$ can be calculated:
90 \frac{d C_{i_\mathrm{ves}}}{dt} = k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]S_\mathrm{ves} -
91 k_{bi}k_{bi\mathrm{adj}}C_{i_\mathrm{ves}}
95 For $k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]$,
96 $k_{fi}$ has units of $\frac{\mathrm{m}}{\mathrm{s}}$,
97 $k_{fi\mathrm{adj}}$ and $k_{bi\mathrm{adj}}$ are unitless,
98 concentration is in units of $\frac{\mathrm{n}}{\mathrm{L}}$, surface
99 area is in units of $\mathrm{m}^2$, $k_{bi}$ has units of
100 $\frac{1}{\mathrm{s}}$ and $C_{i_\mathrm{ves}}$ has units of
101 $\mathrm{n}$, Thus, we have
104 \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} -
105 \frac{1}{\mathrm{s}} \mathrm{n}
107 \frac{\mathrm{m^3}}{\mathrm{s}} \frac{\mathrm{n}}{\mathrm{L}} \frac{1000\mathrm{L}}{\mathrm{m}^3} - \frac{\mathrm{n}}{\mathrm{s}}=
108 \frac{\mathrm{n}}{\mathrm{s}} = 1000 \frac{\mathrm{n}}{\mathrm{s}} - \frac{\mathrm{n}}{\mathrm{s}}
109 \label{eq:state_units}
112 The 1000 isn't in \fref{eq:state} above, because it is unit-dependent.
114 \subsection{Forward adjustments ($k_{fi\mathrm{adj}}$)}
115 \label{sec:forward_adj}
117 The forward rate constant adjustment, $k_{fi\mathrm{adj}}$ takes into
118 account unsaturation ($un_f$), charge ($ch_f$), curvature ($cu_f$),
119 length ($l_f$), and complex formation ($CF1_f$), each of which are
120 modified depending on the specific specie and the vesicle into which
121 the specie is entering.
124 k_{fi\mathrm{adj}} = un_f \cdot ch_f \cdot cu_f \cdot l_f \cdot CF1_f
129 \subsubsection{Unsaturation Forward}
131 In order for a lipid to be inserted into a membrane, a void has to be
132 formed for it to fill. Voids can be generated by the combination of
133 unsaturated and saturated lipids forming herterogeneous domains. Void
134 formation is increased when the unsaturation of lipids in the vesicle
135 is widely distributed; in other words, the insertion of lipids into
136 the membrane is greater when the standard deviation of the
137 unsaturation is larger. Assuming that an increase in width of the
138 distribution linearly decreases the free energy of activation, the
139 $un_f$ parameter must follow
140 $a^{\mathrm{stdev}\left(un_\mathrm{ves}\right)}$ where $a > 1$, so a
141 convenient starting base for $a$ is $2$:
144 un_f = 2^{\mathrm{stdev}\left(un_\mathrm{ves}\right)}
145 \label{eq:unsaturation_forward}
148 The most common $\mathrm{stdev}\left(un_\mathrm{ves}\right)$ is around
149 $1.5$, which leads to a $\Delta \Delta G^\ddagger$ of
150 $\Sexpr{format(digits=3,to.kcal(2^1.5))}
151 \frac{\mathrm{kcal}}{\mathrm{mol}}$.
153 It is not clear that the unsaturation of the inserted monomer will
154 affect the rate of the insertion positively or negatively, so we do
155 not include a term for it in this formalism.
158 \setkeys{Gin}{width=3.2in}
159 <<fig=TRUE,echo=FALSE,results=hide,width=5,height=5>>=
160 curve(2^x,from=0,to=sd(c(0,4)),
161 main="Unsaturation Forward",
162 xlab="Standard Deviation of Unsaturation of Vesicle",
163 ylab="Unsaturation Forward Adjustment")
165 <<fig=TRUE,echo=FALSE,results=hide,width=5,height=5>>=
166 curve(to.kcal(2^x),from=0,to=sd(c(0,4)),
167 main="Unsaturation forward",
168 xlab="Standard Deviation of Unsaturation of Vesicle",
169 ylab="Unsaturation Forward (kcal/mol)")
174 \subsubsection{Charge Forward}
176 A charged lipid such as PS approaching a vesicle with an average
177 charge of the same sign will experience repulsion, whereas those with
178 different signs will experience attraction, the degree of which is
179 dependent upon the charge of the monomer and the average charge of the
180 vesicle. If either the vesicle or the monomer has no charge, there
181 should be no effect of charge upon the rate. This leads us to the
182 following equation, $a^{-\left<ch_v\right> ch_m}$, where
183 $\left<ch_v\right>$ is the average charge of the vesicle, and $ch_m$
184 is the charge of the monomer. If either $\left<ch_v\right>$ or $ch_m$
185 is 0, the adjustment parameter is 1 (no change), whereas it decreases
186 if both are positive or negative, as the product of two real numbers
187 with the same sign is always positive. A convenient base for $a$ is
188 60, resulting in the following equation:
192 ch_f = 60^{-\left<{ch}_v\right> {ch}_m}
193 \label{eq:charge_forward}
196 The most common $\left<{ch}_v\right>$ is around $-0.165$, which leads to
197 a range of $\Delta \Delta G^\ddagger$ from
198 $\Sexpr{format(digits=3,to.kcal(60^(-.165*-1)))}
199 \frac{\mathrm{kcal}}{\mathrm{mol}}$ to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$.
201 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
202 x <- seq(-1,0,length.out=20)
203 y <- seq(-1,0,length.out=20)
204 grid <- expand.grid(x=x,y=y)
205 grid$z <- as.vector(60^(-outer(x,y)))
206 print(wireframe(z~x*y,grid,cuts=50,
208 scales=list(arrows=FALSE),
209 main="Charge Forward",
210 xlab=list("Average Vesicle Charge",rot=30),
211 ylab=list("Component Charge",rot=-35),
212 zlab=list("Charge Forward",rot=93)))
215 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
216 x <- seq(-1,0,length.out=20)
217 y <- seq(-1,0,length.out=20)
218 grid <- expand.grid(x=x,y=y)
219 grid$z <- as.vector(to.kcal(60^(-outer(x,y))))
220 print(wireframe(z~x*y,grid,cuts=50,
222 scales=list(arrows=FALSE),
223 main="Charge Forward (kcal/mol)",
224 xlab=list("Average Vesicle Charge",rot=30),
225 ylab=list("Component Charge",rot=-35),
226 zlab=list("Charge Forward (kcal/mol)",rot=93)))
232 \subsubsection{Curvature Forward}
234 Curvature is a measure of the intrinsic propensity of specific lipids
235 to form micelles (positive curvature), inverted micelles (negative
236 curvature), or planar sheets (zero curvature). In this formalism,
237 curvature is measured as the ratio of the size of the head to that of
238 the base, so negative curvature is bounded by $(0,1)$, zero curvature
239 is 1, and positive curvature is bounded by $(1,\infty)$. The curvature
240 can be transformed into the typical postive/negative mapping using
241 $\log$, which has the additional property of making the range of
242 positive and negative curvature equal, and distributed about 0.
244 As in the case of unsaturation, void formation is increased by the
245 presence of lipids with mismatched curvature. Thus, a larger
246 distribution of curvature in the vesicle increases the rate of lipid
247 insertion into the vesicle. However, a species with curvature $e^{-1}$
248 will cancel out a species with curvature $e$, so we have to log
249 transform (turning these into -1 and 1), then take the absolute value
250 (1 and 1), and finally measure the width of the distribution. Thus, by
251 using the log transform to make the range of the lipid curvature equal
252 between positive and negative, and taking the average to cancel out
253 exactly mismatched curvatures, we come to an equation with the shape
254 $a^{\left<\log cu_\mathrm{vesicle}\right>}$. A convenient base for $a$
259 % cu_f = 10^{\mathrm{stdev}\left|\log cu_\mathrm{vesicle}\right|}
260 cu_f = 10^{\left|\left<\log cu_\mathrm{vesicle} \right>\right|\mathrm{stdev} \left|\log cu_\mathrm{vesicle}\right|}
261 \label{eq:curvature_forward}
264 The most common $\left|\left<\log {cu}_v\right>\right|$ is around
265 $0.013$, which with the most common $\mathrm{stdev} \log
266 cu_\mathrm{vesicle}$ of $0.213$ leads to a $\Delta \Delta G^\ddagger$
267 of $\Sexpr{format(digits=3,to.kcal(10^(0.13*0.213)))}
268 \frac{\mathrm{kcal}}{\mathrm{mol}}$. This is a consequence of the
269 relatively matched curvatures in our environment.
271 % 1.5 to 0.75 3 to 0.33
272 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
273 grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
274 y sd(abs(log(c(1,0.33)))),sd(abs(log(c(0.33,3)))))),length.out=20),
275 y=seq(0,max(c(mean(log(c(1,3)),
276 mean(log(c(1,0.33))),
277 mean(log(c(0.33,3)))))),length.out=20))
278 grid$z <- 10^(grid$x*grid$y)
279 print(wireframe(z~x*y,grid,cuts=50,
281 scales=list(arrows=FALSE),
282 xlab=list("Vesicle stdev log curvature",rot=30),
283 ylab=list("Vesicle average log curvature",rot=-35),
284 zlab=list("Vesicle Curvature Forward",rot=93)))
287 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
288 grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
289 sd(abs(log(c(1,0.33)))),sd(abs(log(c(0.33,3)))))),length.out=20),
290 y=seq(0,max(c(mean(log(c(1,3)),
291 mean(log(c(1,0.33))),
292 mean(log(c(0.33,3)))))),length.out=20))
293 grid$z <- to.kcal(10^(grid$x*grid$y))
294 print(wireframe(z~x*y,grid,cuts=50,
296 scales=list(arrows=FALSE),
297 xlab=list("Vesicle stdev log curvature",rot=30),
298 ylab=list("Vesicle average log curvature",rot=-35),
299 zlab=list("Vesicle Curvature Forward (kcal/mol)",rot=93)))
304 \subsubsection{Length Forward}
306 As in the case of unsaturation, void formation is easier when vesicles
307 are made up of components of widely different lengths. Thus, when the
308 width of the distribution of lengths is larger, the forward rate
309 should be greater as well, leading us to an equation of the form
310 $x^{\mathrm{stdev} l_\mathrm{ves}}$, where $\mathrm{stdev}
311 l_\mathrm{ves}$ is the standard deviation of the length of the
312 components of the vesicle, which has a maximum possible value of 8 and
313 a minimum of 0 in this set of experiments. A convenient base for $x$
317 l_f = 2^{\mathrm{stdev} l_\mathrm{ves}}
318 \label{eq:length_forward}
321 The most common $\mathrm{stdev} l_\mathrm{ves}$ is around $3.4$, which leads to
322 a range of $\Delta \Delta G^\ddagger$ of
323 $\Sexpr{format(digits=3,to.kcal(2^(3.4)))}
324 \frac{\mathrm{kcal}}{\mathrm{mol}}$.
326 While it could be argued that increased length of the monomer could
327 affect the rate of insertion into the membrane, it's not clear whether
328 it would increase (by decreasing the number of available hydrogen
329 bonds, for example) or decrease (by increasing the time taken to fully
330 insert the acyl chain, for example) the rate of insertion or to what
331 degree, so we do not take it into account in this formalism.
333 \fixme{Incorporate McLean84 here}
334 From McLean84LIB: Although it is difficult to measure cmc values for
335 the sparingly soluble lipids used in this study, estimates for
336 lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds,
337 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution;
338 Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6
339 X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X
340 10-8 M for DMPC was estimated from the linear relationship between ln
341 cmc and the number of carbons in the PC acyl chain by using data for n
342 = 7, 8, 10, and 16 [summarized in Tanford (1980)].
344 From Nichols85: The magnitude of the dissociation rate constant
345 decreases by a factor of approximately 3.2 per carbon increase in acyl
346 chain length (see Table II here) {take into account for the formula;
349 From Nichols85: The magnitude of the partition coefficient increases
350 with acyl chain length [Keq(M-C6-NBD-PC) = (9.8±2.1) X 106 M-1 and Keq
351 (P-C6-NBD-PC) = (9.4±1.3) x 107 M-1. {take into account for the
352 formula; rz 8/17/2010}.
354 From Nichols85: The association rate constant is independent of acyl
355 chain length. {take into account for the formula; rz 8/17/2010}.
358 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=5>>=
359 curve(2^x,from=0,to=sd(c(12,24)),
360 main="Length forward",
361 xlab="Standard Deviation of Length of Vesicle",
362 ylab="Length Forward Adjustment")
364 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=5>>=
365 curve(to.kcal(2^x),from=0,to=sd(c(12,24)),
366 main="Length forward",
367 xlab="Standard Deviation of Length of Vesicle",
368 ylab="Length Forward Adjustment (kcal/mol)")
372 \subsubsection{Complex Formation}
373 There is no contribution of complex formation to the forward reaction
374 rate in the current formalism.
378 \label{eq:complex_formation_forward}
381 \subsection{Backward adjustments ($k_{bi\mathrm{adj}}$)}
383 Just as the forward rate constant adjustment $k_{fi\mathrm{adj}}$
384 does, the backwards rate constant adjustment $k_{bi\mathrm{adj}}$
385 takes into account unsaturation ($un_b$), charge ($ch_b$), curvature
386 ($cu_b$), length ($l_b$), and complex formation ($CF1_b$), each of
387 which are modified depending on the specific specie and the vesicle
388 into which the specie is entering:
392 k_{bi\mathrm{adj}} = un_b \cdot ch_b \cdot cu_b \cdot l_b \cdot CF1_b
396 \subsubsection{Unsaturation Backward}
398 Unsaturation also influences the ability of a lipid molecule to leave
399 a membrane. If a molecule has an unsaturation level which is different
400 from the surrounding membrane, it will be more likely to leave the
401 membrane. The more different the unsaturation level is, the greater
402 the propensity for the lipid molecule to leave. However, a vesicle
403 with some unsaturation is more favorable for lipids with more
404 unsaturation than the equivalent amount of less unsatuturation, so the
405 difference in energy between unsaturation is not linear. Therefore, an
406 equation with the shape
407 $x^{\left| y^{-\left< un_\mathrm{ves}\right> }-y^{-un_\mathrm{monomer}} \right| }$
408 where $\left<un_\mathrm{ves}\right>$ is the average unsaturation of
409 the vesicle, and $un_\mathrm{monomer}$ is the average unsaturation. In
410 this equation, as the average unsaturation of the vesicle is larger,
413 un_b = 7^{1-\left(20\left(2^{-\left<un_\mathrm{vesicle} \right>} - 2^{-un_\mathrm{monomer}} \right)^2+1\right)^{-1}}
414 \label{eq:unsaturation_backward}
417 The most common $\left<un_\mathrm{ves}\right>$ is around $1.7$, which leads to
418 a range of $\Delta \Delta G^\ddagger$ from
419 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-0)^2+1))))}
420 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with 0 unsaturation
422 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-4)^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
423 for monomers with 4 unsaturations.
426 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
427 grid <- expand.grid(x=seq(0,4,length.out=20),
428 y=seq(0,4,length.out=20))
429 grid$z <- (7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1)))
430 print(wireframe(z~x*y,grid,cuts=50,
432 scales=list(arrows=FALSE),
433 xlab=list("Average Vesicle Unsaturation",rot=30),
434 ylab=list("Monomer Unsaturation",rot=-35),
435 zlab=list("Unsaturation Backward",rot=93)))
438 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
439 grid <- expand.grid(x=seq(0,4,length.out=20),
440 y=seq(0,4,length.out=20))
441 grid$z <- to.kcal((7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1))))
442 print(wireframe(z~x*y,grid,cuts=50,
444 scales=list(arrows=FALSE),
445 xlab=list("Average Vesicle Unsaturation",rot=30),
446 ylab=list("Monomer Unsaturation",rot=-35),
447 zlab=list("Unsaturation Backward (kcal/mol)",rot=93)))
454 \subsubsection{Charge Backwards}
455 As in the case of monomers entering a vesicle, monomers leaving a
456 vesicle leave faster if their charge has the same sign as the average
457 charge vesicle. An equation of the form $ch_b = a^{\left<ch_v\right>
458 ch_m}$ is then appropriate, and using a base of $a=20$ yields:
461 ch_b = 20^{\left<{ch}_v\right> {ch}_m}
462 \label{eq:charge_backwards}
465 The most common $\left<ch_v\right>$ is around $-0.164$, which leads to
466 a range of $\Delta \Delta G^\ddagger$ from
467 $\Sexpr{format(digits=3,to.kcal(20^(-.164*-1)))}
468 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with charge $-1$ to
469 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$
470 for monomers with charge $0$.
473 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
474 x <- seq(-1,0,length.out=20)
475 y <- seq(-1,0,length.out=20)
476 grid <- expand.grid(x=x,y=y)
477 grid$z <- as.vector(20^(outer(x,y)))
478 print(wireframe(z~x*y,grid,cuts=50,
480 scales=list(arrows=FALSE),
481 xlab=list("Average Vesicle Charge",rot=30),
482 ylab=list("Component Charge",rot=-35),
483 zlab=list("Charge Backwards",rot=93)))
486 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
487 x <- seq(-1,0,length.out=20)
488 y <- seq(-1,0,length.out=20)
489 grid <- expand.grid(x=x,y=y)
490 grid$z <- to.kcal(as.vector(20^(outer(x,y))))
491 print(wireframe(z~x*y,grid,cuts=50,
493 scales=list(arrows=FALSE),
494 xlab=list("Average Vesicle Charge",rot=30),
495 ylab=list("Component Charge",rot=-35),
496 zlab=list("Charge Backwards (kcal/mol)",rot=93)))
501 \subsubsection{Curvature Backwards}
503 The less a monomer's intrinsic curvature matches the average curvature
504 of the vesicle in which it is in, the greater its rate of efflux. If
505 the difference is 0, $cu_f$ needs to be one. To map negative and
506 positive curvature to the same range, we also need take the logarithm.
507 Increasing mismatches in curvature increase the rate of efflux, but
508 asymptotically. An equation which satisfies this critera has the
509 form $cu_f = a^{1-\left(b\left( \left< \log cu_\mathrm{vesicle} \right>
510 -\log cu_\mathrm{monomer}\right)^2+1\right)^{-1}}$. An
511 alternative form would use the aboslute value of the difference,
512 however, this yields a cusp and sharp increase about the curvature
513 equilibrium, which is decidedly non-elegant. We have chosen bases of
517 cu_b = 7^{1-\left(20\left(\left< \log cu_\mathrm{vesicle} \right> -\log cu_\mathrm{monomer} \right)^2+1\right)^{-1}}
518 \label{eq:curvature_backwards}
521 The most common $\left<\log cu_\mathrm{ves}\right>$ is around $-0.013$, which leads to
522 a range of $\Delta \Delta G^\ddagger$ from
523 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(0.8))^2+1))))}
524 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature 0.8
526 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(1.3))^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
527 for monomers with curvature 1.3 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature near 1.
530 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
531 grid <- expand.grid(x=seq(0.8,1.33,length.out=20),
532 y=seq(0.8,1.33,length.out=20))
533 grid$z <- 7^(1-1/(20*(log(grid$x)-log(grid$y))^2+1))
534 print(wireframe(z~x*y,grid,cuts=50,
536 scales=list(arrows=FALSE),
537 xlab=list("Vesicle Curvature",rot=30),
538 ylab=list("Monomer Curvature",rot=-35),
539 zlab=list("Curvature Backward",rot=93)))
542 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
543 grid <- expand.grid(x=seq(0.8,1.33,length.out=20),
544 y=seq(0.8,1.33,length.out=20))
545 grid$z <- to.kcal(7^(1-1/(20*(log(grid$x)-log(grid$y))^2+1)))
546 print(wireframe(z~x*y,grid,cuts=50,
548 scales=list(arrows=FALSE),
549 xlab=list("Vesicle Curvature",rot=30),
550 ylab=list("Monomer Curvature",rot=-35),
551 zlab=list("Curvature Backward (kcal/mol)",rot=93)))
556 \subsubsection{Length Backwards}
558 In a model membrane, the dissociation constant increases by a factor
559 of approximately 3.2 per carbon decrease in acyl chain length (Nichols
560 1985). Unfortunatly, the experimental data known to us only measures
561 chain length less than or equal to the bulk lipid, and does not exceed
562 it, and is only known for one bulk lipid species (DOPC). We assume
563 then, that the increase is in relationship to the average vesicle, and
564 that lipids with larger acyl chain length will also show an increase
565 in their dissociation constant.
568 l_b = 3.2^{\left|\left<l_\mathrm{ves}\right>-l_\mathrm{monomer}\right|}
569 \label{eq:length_backward}
572 The most common $\left<\log l_\mathrm{ves}\right>$ is around $17.75$, which leads to
573 a range of $\Delta \Delta G^\ddagger$ from
574 $\Sexpr{format(digits=3,to.kcal(3.2^abs(12-17.75)))}
575 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with length 12
577 $\Sexpr{format(digits=3,to.kcal(3.2^abs(24-17.75)))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
578 for monomers with length 24 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature near 18.
581 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
582 grid <- expand.grid(x=seq(12,24,length.out=20),
583 y=seq(12,24,length.out=20))
584 grid$z <- 3.2^(abs(grid$x-grid$y))
585 print(wireframe(z~x*y,grid,cuts=50,
587 scales=list(arrows=FALSE),
588 xlab=list("Average Vesicle Length",rot=30),
589 ylab=list("Monomer Length",rot=-35),
590 zlab=list("Length Backward",rot=93)))
593 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
594 grid <- expand.grid(x=seq(12,24,length.out=20),
595 y=seq(12,24,length.out=20))
596 grid$z <- to.kcal(3.2^(abs(grid$x-grid$y)))
597 print(wireframe(z~x*y,grid,cuts=50,
599 scales=list(arrows=FALSE),
600 xlab=list("Average Vesicle Length",rot=30),
601 ylab=list("Monomer Length",rot=-35),
602 zlab=list("Length Backward (kcal/mol)",rot=93)))
608 \subsubsection{Complex Formation Backward}
610 Complex formation describes the interaction between CHOL and PC or SM,
611 where PC or SM protects the hydroxyl group of CHOL from interactions
612 with water, the ``Umbrella Model''. PC ($CF1=2$) can interact with two
613 CHOL, and SM ($CF1=3$) with three CHOL ($CF1=-1$). If the average of
614 $CF1$ is positive (excess of SM and PC with regards to complex
615 formation), species with negative $CF1$ (CHOL) will be retained. If
616 average $CF1$ is negative, species with positive $CF1$ are retained.
617 An equation which has this property is
618 $CF1_b=a^{\left<CF1_\mathrm{ves}\right>
619 CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{ves}\right>
620 CF1_\mathrm{monomer}\right|}$, where difference of the exponent is
621 zero if the average $CF1$ and the $CF1$ of the specie have the same
622 sign, or double the product if the signs are different. A convenient
623 base for $a$ is $1.5$.
627 CF1_b=1.5^{\left<CF1_\mathrm{ves}\right> CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{ves}\right> CF1_\mathrm{monomer}\right|}
628 \label{eq:complex_formation_backward}
631 The most common $\left<CF1_\mathrm{ves}\right>$ is around $0.925$,
632 which leads to a range of $\Delta \Delta G^\ddagger$ from
633 $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*-1-abs(0.925*-1))))}
634 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
636 $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*2-abs(0.925*2))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
637 for monomers with complex formation $2$ to
638 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
642 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
643 grid <- expand.grid(x=seq(-1,3,length.out=20),
644 y=seq(-1,3,length.out=20))
645 grid$z <- 1.5^(grid$x*grid$y-abs(grid$x*grid$y))
646 print(wireframe(z~x*y,grid,cuts=50,
648 scales=list(arrows=FALSE),
649 xlab=list("Vesicle Complex Formation",rot=30),
650 ylab=list("Monomer Complex Formation",rot=-35),
651 zlab=list("Complex Formation Backward",rot=93)))
654 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
655 grid <- expand.grid(x=seq(-1,3,length.out=20),
656 y=seq(-1,3,length.out=20))
657 grid$z <- to.kcal(1.5^(grid$x*grid$y-abs(grid$x*grid$y)))
658 print(wireframe(z~x*y,grid,cuts=50,
660 scales=list(arrows=FALSE),
661 xlab=list("Vesicle Complex Formation",rot=30),
662 ylab=list("Monomer Complex Formation",rot=-35),
663 zlab=list("Complex Formation Backward (kcal/mol)",rot=93)))
669 \subsection{Per-Lipid Kinetic Parameters}
671 Each of the 5 lipid types have different kinetic parameters; to the
672 greatest extent possible, we have derived these from literature.
676 \begin{tabular}{c c c c c c c}
677 Type & $k_f$ & $k_b$ & Area (\r{A}$^2$) & Charge & CF1 & Curvature \\
679 PC & $3.7\cdot 10^6$ & $2\cdot 10^{-5}$ & 63 & 0 & 2 & 0.8 \\
680 PS & $3.7\cdot 10^6$ & $1.5\cdot 10^{-5}$ & 54 & -1 & 0 & 1 \\
681 CHOL & $5.1\cdot 10^7$ & $2.8\cdot 10^{-4}$ & 38 & 0 & -1 & 1.21 \\
682 SM & $3.7\cdot 10^6$ & $3.1\cdot 10^{-3}$ & 51 & 0 & 3 & 0.8 \\
683 PE & $2.3\cdot 10^6$ & $10^{-5}$ & 55 & 0 & 0 & 1.33 \\
685 \caption{Kinetic parameters of lipid types}
686 \label{tab:kinetic_parameters_lipid_types}
689 \subsubsection{$k_f$ for lipid types}
690 For PC, $k_f$ was measured by Nichols85 to be $3.7\cdot 10^6
691 \frac{1}{\mathrm{M}\cdot \mathrm{s}}$ by the partitioning of
692 P-C$_6$-NBD-PC between DOPC vesicles and water. The method utilized by
693 Nichols85 has the weakness of using NBD-PC, with associated label
694 perturbations. As similar measures do not exist for SM or PS, we
695 assume that they have the same $k_f$. For CHOL, Estronca07 found a
696 value for $k_f$ of $5.1\cdot 10^7 \frac{1}{\mathrm{M}\cdot
697 \mathrm{s}}$. For PE, Abreu04 found a value for $k_f$ of $2.3\cdot
698 10^6$. \fixme{I'm missing the notes on these last two papers, so this
701 \subsubsection{$k_b$ for lipid types}
703 $k_b$ for PC was measured by Wimley90 using a radioactive label and
704 large unilammelar vesicles at 30\textdegree C. The other values were
705 calculated from the experiments of Nichols82 where the ratio of $k_b$
706 of different types was measured to that of PC.
707 See~\fref{tab:kinetic_parameters_lipid_types}.
709 assigned accordingly. kb(PS) was assumed to be the same as kb(PG)
710 given by Nichols82 (also ratio from kb(PC)). kb(SM) is taken from
711 kb(PC) of Wimley90 (radioactive), and then a ratio of kb(PC)/kb(SM)
712 taken from Bai97: = 34/2.2 = 15.45; 2.0 x 10-4 x 15.45 = 3.1 x 10-3 s
713 -1. kb(CHOL) taken from Jones90 (radioactive; POPC LUV; 37°).
716 PE 0.45 <- from Nichols82
722 \subsubsection{Area for lipid types}
725 From Sampaio05: Besides this work and our own earlier report on the
726 association of NBD-DMPE with lipid bilayers (Abreu et al., 2004), we
727 are aware of only one other report in the literature (Nichols, 1985)
728 in which all the kinetic constants of lipid-derived amphiphile
729 association with lipid bilayer membranes were experimentally measured.
730 {indeed; everything is k- !!!; rz}
732 From McLean84LIB: Although it is difficult to measure cmc values for
733 the sparingly soluble lipids used in this study, estimates for
734 lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds,
735 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution;
736 Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6
737 X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X
738 10-8 M for DMPC was estimated from the linear relationship between ln
739 cmc and the number of carbons in the PC acyl chain by using data for n
740 = 7, 8, 10, and 16 [summarized in Tanford (1980)].
742 From Toyota08: Recently, several research groups have reported
743 self-reproducing systems of giant vesicles that undergo a series of
744 sequential transformations: autonomous growth, self-division, and
745 chemical reactions to produce membrane constituents within the giant
748 Vesicle sizes of 25 nm for SUV and 150 nm for LUV were mentioned by
751 From Lund-Katz88: Charged and neutral small unilamellar vesicles
752 composed of either saturated PC, unsaturated PC, or SM had similar
753 size distributions with diameters of 23 \& 2 nm.
755 From Sampaio05LIB: The exchange of lipids and lipid derivatives
756 between lipid bilayer vesicles has been studied for at least the last
757 30 years. Most of this work has examined the exchange of amphiphilic
758 molecules between a donor and an acceptor population. The measured
759 efflux rates were shown in almost all cases, not surprisingly, to be
760 first order processes. In all of this work, the insertion rate has
761 been assumed to be much faster than the efflux rate. Having measured
762 both the insertion and desorption rate constants for amphiphile
763 association with membranes, our results show that this assumption is
764 valid. In several cases reported in the literature, the insertion rate
765 constant was assumed, although never demonstrated, to be a
766 diffusion-controlled process.
768 (for methods? From McLean84LIB: The activation free energies and free
769 energies of transfer from self-micelles to water increase by 2.2 and
770 2.1 kJ mol-' per methylene group, respectively. {see if we can use it
771 to justify arranging our changed in activating energy around 1
774 Jones90 give diameter of LUV as 100 nm, and of SUV as 20 nm; that
775 would give the number of molecules per outer leaflet of a vesicle as
778 Form Simard08: Parallel studies with SUV and LUV revealed that
779 although membrane curvature does have a small effect on the absolute
780 rates of FA transfer between vesicles, the ΔG of membrane desorption
781 is unchanged, suggesting that the physical chemical properties which
782 govern FA desorption are dependent on the dissociating molecule rather
783 than on membrane curvature. However, disagreements on this fundamental
784 issue continue (57, 61, 63, 64)
786 (methods regarding the curvature effect: Kleinfeld93 showed that the
787 transfer parameters of long-chain FFA between the lipid vesicles
788 depend on vesicle curvature and composition. Transfer of stearic acid
789 is much slower from LUV as compared to SUV).
791 From McLean84: In a well-defined experimental system consisting of
792 unilamellar lipid vesicles, in the absence of protein, the
793 rate-limiting step for the overall exchange process is desorption
794 (McLean \& Phillips, 1981). {thus I can take exchange data for the
795 estimation of k- rz; 8/11/08}.
797 \subsubsection{Complex Formation 1}
799 From Thomas88a: SM decreases the rate of cholesterol transfer, while
800 phosphatidylethanolamine (PE) and phosphatidylserine (PS) have no
801 effect at physiologically significant levels.
804 \section{Simulation Methodology}
806 \subsection{Overall Architecture}
808 The simulation is currently run by single program written in perl
809 using various different modules to handle the subsidiary parts. It
810 produces output for each generation, including the step immediately
811 preceeding and immediately following a vesicle split (and optionally,
812 each step) that is written to a state file which contains the state of
813 the vesicle, environment, kinetic parameters, program invocation
814 options, time, and various other parameters necessary to recreate the
815 state vector at a given time. This output file is then read by a
816 separate program in perl to produce different output which is
817 generated out-of-band; output which includes graphs and statistical
818 analysis is performed using R (and various grid graphics modules)
819 called from the perl program.
821 The separation of simulation and output generation allows refining
822 output, and simulations performed on different versions of the
823 underlying code to be compared using the same analysis methods and
824 code. It also allows later simulations to be restarted from a specific
825 generation, utilizing the same environment.
827 Random variables of different distributions are calculated using
828 Math::Random, which is seeded for each run using entropy from the
829 linux kernel's urandom device.
831 \subsection{Environment Creation}
833 \subsubsection{Components}
834 The environment contains different concentrations of different
835 components. In the current set of simulations, there are
836 \Sexpr{1+4*5*7} different components, consisting of PC, PE, PS, SM,
837 and CHOL, with all lipids except for CHOL having 5 possible
838 unsaturations rangiong from 0 to 4, and 7 lengths from $12,14,...,22$
839 ($7\cdot 5\cdot4+1=\Sexpr{1+4*5*7}$). In cases where the environment
840 has less than the maximum number of components, the components are
841 selected in order without replacement from a randomly shuffled deck of
842 component (with the exception of CHOL) represented once until the
843 desired number of unique components are obtained. CHOL is over
844 representated $\Sexpr{5*7}$ times to be at the level of other lipid
845 types, ensures that the probability of CHOL being asbent in the
846 environment is the same as the probability of one of the other lipid
847 types (PS, PE, etc.) being entirely absent. This reduces the number of
848 simulations with a small number of components which are entirely
851 \subsubsection{Concentration}
852 Once the components of the environment have been selected, the
853 concentration of thoes components is determined. In experiments where
854 the environmental concentration is the same across all lipid
855 components, the concentration is set to $10^{-10}\mathrm{M}$. In other
856 cases, the environmental concentration is set to a random number from
857 a gamma distribution with shape parameter 2 and an average of
858 $10^{-10}\mathrm{M}$. The base concentration ($10^{-10}\mathrm{M}$)
859 can also be altered in the initialization of the experiment to
860 specific values for specific lipid types or components.
862 \subsection{Initial Vesicle Creation}
864 Initial vesicles are seeded by selecting lipid molecules from the
865 environment until the vesicle reaches a specific starting size. The
866 vesicle starting size has gamma distribution with shape parameter 2
867 and a mean of the experimentally specified starting size, with a
868 minimum of 5 lipid molecules. Lipid molecules are then selected to be
869 added to the lipid membrane according to four different methods. In
870 the constant method, molecules are added in direct proportion to their
871 concentration in the environment. The uniform method adds molecules in
872 proportion to their concentration in the environment scaled by a
873 uniform random value, whereas the random method adds molecules in
874 proportion to a uniform random value. The final method is a binomial
875 method, which adjust the porbability of adding a molecule of a
876 specific component by the concentration of that component, and then
877 adds the molecules one by one to the membrane. This final method is
878 also used in the first three methods to add any missing molecules to
879 the starting vesicle which were unallocated due to the requirement to
880 add integer numbers of molecules.
882 \subsection{Simulation Step}
884 Once the environment has been created and the initial vesicle has been
885 formed, molecules join and leave the vesicle based on the kinetic
886 parameters and state equation discussed until the vesicle splits
887 forming two daughter vesicles, one of which the program continues to
890 \subsubsection{Calculation of Vesicle Properties}
893 $S_\mathrm{ves}$ is the surface area of the vesicle, and is the sum of
894 the surface area of all of the individual lipid components; each lipid
895 type has a different surface area; we na\"ively assume that the lipid
896 packing is optimal, and there is no empty space.
898 \subsubsection{Joining and Leaving of Lipid Molecules}
900 Determining the number of molecules to add to the lipid membrane
901 ($n_i$) requires knowing $k_{fi\mathrm{adj}}$, the surface area of the
902 vesicle $S_\mathrm{ves}$ (see~\fref{sec:ves_prop}), the time interval
903 $dt$ during which lipids are added, the base $k_{fi}$, and the
904 concentration of the monomer in the environment
905 $\left[C_{i\mathrm{ves}}\right]$ (see~\fref{eq:state}).
906 $k_{fi\mathrm{adj}}$ is calculated (see~\fref{eq:kf_adj}) based on the
907 vesicle properties and their hypothesized effect on the rate (in as
908 many cases as possible, experimentally based)
909 (see~\fref{sec:forward_adj} for details). $dt$ can be varied
910 (see~\fref{sec:step_duration}), but for a given step is constant. This
911 leads to the following:
913 $n_i = k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]S_\mathrm{ves}\mathrm{N_A}dt$
915 In the cases where $n_i > 1$, the integer number of molecules is
916 added. Fractional $n_i$ or the fractional remainder after the addition
917 of the integer molecules are the probability of adding a specific
918 molecule, and are compared to a uniformly distributed random value
919 between 0 and 1. If the random value is less than or equal to the
920 fractional part of $n_i$, an additional molecule is added.
922 Molecules leaving the vesicle are handled in a similar manner, with
924 $n_i = k_{bi}k_{bi\mathrm{adj}}C_{i_\mathrm{ves}}\mathrm{N_A}dt$.
926 While programatically, the molecule removal happens after the
927 addition, the properties that each operates on are the same, so they
928 can be considered to have been added and removed at the same instant.
929 [This also avoids cases where a removal would have resulted in
930 negative molecules for a particular type, which is obviously
933 \subsubsection{Step duration}
934 \label{sec:step_duration}
936 $dt$ is the time taken for each step of joining, leaving, and checking
937 split. In the current implementation, it starts with a value of
938 $10^{-6}\mathrm{s}$ but this is modified in between each step if the
939 number of molecules joining or leaving is too large or small. If more
940 than half of the starting vesicle size molecules join or leave in a
941 single step, $dt$ is reduced by half. If less than the starting
942 vesicle size molecules join or leave in 100 steps, $dt$ is doubled.
943 This is necessary to curtail run times and to automatically adjust to
944 different experimental runs.
946 (In every run seen so far, the initial $dt$ was too small, and was
947 increased before the first generation occured; at no time was $dt$ too
950 \subsubsection{Vesicle splitting}
952 If a vesicle has grown to a size which is more than double the
953 starting vesicle size, the vesicle splits. More elaborate mechanisms
954 for determining whether a vesicle should split are of course possible,
955 but not currently implemented. Molecules are assigned to the daughter
956 vesicles at random, with each daughter vesicle having an equal chance
957 of getting a single molecule. The number of molecules to assign to the
958 first vesicle has binomial distribution with a probability of an event
959 in each trial of 0.5 and a number of trials equal to the number of
964 The environment, initial vesicle, and the state of the vesicle
965 immediately before and immediately after splitting are stored to disk
966 to produce later output.
968 \section{Analyzing output}
970 Analyzing of output is handled by a separate perl program which shares
971 many common modules with the simulation program. Current output
972 includes simulation progress, summary tables, summary statistics, and
975 \subsection{PCA plots}
977 Vesicles have many different axes which contribute to their variation
978 between subsequent generations; two major groups of axes are the
979 components and properties of vesicles. Each component in a vesicle is
980 an axis on its own; it can be measured either as an absolute number of
981 molecules in each component, or the fraction of molecules of that
982 component over the total number of molecules; the second approach is
983 often more convenient, as it allows vesicles of different number of
984 molecules to be more directly compared (though it hides any affect of
987 In order to visualize the transition of subsequent generations of
988 vesicles from their initial state in the simulation, to their final
989 state at the termination of
991 \subsection{Carpet plots}
995 % \bibliographystyle{plainnat}
996 % \bibliography{references.bib}