1 \documentclass[english,12pt]{article}
3 \setmainfont{FreeSerif}
10 %\usepackage[pdftex]{graphicx}
12 \usepackage[bf]{caption}
15 %\usepackage{textcomp}
16 %\usepackage{mathrsfs}
20 %\usepackage[light,all]{draftcopy}
25 \usepackage[backend=biber,natbib=true,hyperref=true,citestyle=numeric-comp,style=nature,autocite=inline]{biblatex}
26 \addbibresource{references.bib}
27 \usepackage[hyperfigures,bookmarks,colorlinks,citecolor=black,filecolor=black,linkcolor=black,urlcolor=black]{hyperref}
28 \usepackage[noblocks,auth-sc]{authblk}
29 \usepackage[capitalise]{cleveref}
30 \usepackage[markifdraft,raisemark=0.01\paperheight]{gitinfo2}
31 %\usepackage[sectionbib,sort&compress,numbers]{natbib}
32 % \usepackage[nomargin,inline,draft]{fixme}
33 %\usepackage[x11names,svgnames]{xcolor}
34 %\usepackage{texshade}
35 \renewcommand{\textfraction}{0.15}
36 \renewcommand{\topfraction}{0.85}
37 \renewcommand{\bottomfraction}{0.65}
38 \renewcommand{\floatpagefraction}{0.60}
39 %\renewcommand{\baselinestretch}{1.8}
41 \newcommand{\DLA}[1]{\textcolor{red}{\fxnote{DLA: #1}}}
43 \newcommand{\RZ}[1]{\textcolor{red}{\fxnote{RZ: #1}}}
45 \newcommand{\OM}[1]{\textcolor{red}{\fxnote{OM: #1}}}
55 \title{Appendix 1: Kinetic formalism for R-GARD Simulations}
57 % \author[1,2]{Don Armstrong}
58 % %\ead{don@donarmstrong.com}
59 % \author[2]{Raphael Zidovetzki}
60 % %\ead{raphael.zidovetzki@ucr.edu}
61 % \affil[1]{Carl R. Woese Institute for Genomic Biology, University of Illinois at Urbana-Champaign, Urbana, IL, USA}
62 % \affil[2]{Cell Biology and Neuroscience, University of California at Riverside, Riverside, CA, USA}
63 % \affil[ ]{don@donarmstrong.com, raphael.zidovetzki@ucr.edu}
69 <<load.libraries,echo=FALSE,results="hide",warning=FALSE,message=FALSE,error=FALSE,cache=FALSE>>=
70 opts_chunk$set(dev="CairoPDF",out.width="\\columnwidth",out.height="0.7\\textheight",out.extra="keepaspectratio")
71 opts_chunk$set(cache=TRUE, autodep=TRUE)
72 options(scipen = -1, digits = 2)
77 to.latex <- function(x){
78 gsub("\\\\","\\\\\\\\",latexSN(x))
81 to.kcal <- function(k,temp=300) {
83 return(-gasconst*temp*log(k)/1000)
87 % \setcounter{figure}{0} \setcounter{table}{0}
89 \renewcommand{\thefigure}{S\@arabic\c@figure}
90 \renewcommand{\thetable}{S\@arabic\c@table}
93 % \section{Competition Implementation}
94 % \subsection{Implementation changes}
97 % \item settable maximum number of vesicles to track (default $10^4$)
98 % \item start with 1~L ($10^{-3}$~m$^3$) cube
99 % \item if at any point the number of vesicles exceeds the maximum
100 % number, chop the volume and environment molecule number into tenths,
101 % randomly select one tenth of the vesicles, and continue tracking.
102 % \item generations will be counted per vesicle, and each progeny
103 % vesicle will have a generation number one greater than the parental
105 % \item 100 generations can result in as many as $2^{100}$
106 % ($\Sexpr{2^100}$) vesicles or as few as
108 % \item Environment will use a specific number of each component instead
109 % of a constant concentration; as the number may be larger than
110 % \texttt{long long} ($2^{64}$), we use libgmp to handle an arbitrary
111 % precision number of components
114 % \subsection{Infrastructure changes}
116 % \item Rewrite core bits in C
117 % \item Use libgmp for handling large ints
118 % \item Use openmpi to split the calculations out over multiple
119 % machines/processors and allow deploying to larger
120 % clusters/supercomputers
125 \section{State Equation}
126 % double check this with the bits in the paper
128 For a system with monomers $(_\mathrm{monomer})$ and a vesicle
129 $(_\mathrm{vesicle})$, the change in concentration of the $i$th component of
130 a lipid vesicle per change in time ($d \left[C_{i_\mathrm{vesicle}}\right]/dt$)
131 can be described by a modification of the basic mass action law:
134 \frac{d \left[C_{i_\mathrm{vesicle}}\right]}{dt} = k_{\mathrm{f}i}k_{\mathrm{f}i\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]\left[S_\mathrm{vesicle}\right] -
135 k_{\mathrm{b}i}k_{\mathrm{b}i\mathrm{adj}}\left[C_{i_\mathrm{vesicle}}\right]
139 The base forward kinetic parameter for the $i$th component is $k_{\mathrm{f}i}$
140 and is dependent on the particular lipid type (\ac{PC}, \ac{PS}, \ac{SM}, etc.). The
141 forward adjustment parameter, $k_{\mathrm{f}i\mathrm{adj}}$, is based on the
142 properties of the vesicle and the specific component (type, length,
143 unsaturation, etc.) (see \cref{eq:kf_adj,sec:kinetic_adjustments}).
144 $\left[C_{i_\mathrm{monomer}}\right]$ is the molar concentration of
145 monomer of the $i$th component. $\left[S_\mathrm{vesicle}\right]$ is the surface
146 area of the vesicle per volume. The base backwards kinetic parameter
147 for the $i$th component is $k_{\mathrm{b}i}$ and its adjustment parameter
148 $k_{\mathrm{b}i\mathrm{adj}}$ (see \cref{eq:kb_adj,sec:kinetic_adjustments}).
149 $\left[C_{i_\mathrm{vesicle}}\right]$ is the molar concentration of
150 the $i$th component in the vesicle.
152 \subsection{Per-Lipid Kinetic Parameters}
154 <<kf_prime,echo=FALSE,results='hide'>>=
155 kf.prime <- c(3.7e6,3.7e6,5.1e7,3.7e6,2.3e6)
156 kf <- (as.numeric(kf.prime)*10^-3)/(63e-20*6.022e23)
159 Each of the 5 lipid types has different kinetic parameters; where
160 available, these were taken from literature (\cref{tab:kinetic_parameters_lipid_types}).
164 \begin{tabular}{c c c c c c c c}
166 Type & $k_\mathrm{f}$ $\left(\frac{\mathrm{m}}{\mathrm{s}}\right)$
167 & $k'_\mathrm{f}$ $\left(\frac{1}{\mathrm{M} \mathrm{s}}\right)$
168 & $k_\mathrm{b}$ $\left(\mathrm{s}^{-1}\right)$
169 & Area $\left(\mathrm{Å}^2\right)$ & Charge & $\mathrm{CF}1$ & Curvature \\
171 PC & $\Sexpr{kf[1]}$ & $3.7 \times 10^6$ & $2 \times 10^{-5}$ & 63 & 0 & 2 & 0.8 \\
172 PS & $\Sexpr{kf[2]}$ & $3.7 \times 10^6$ & $1.25\times 10^{-5}$ & 54 & -1 & 0 & 1 \\
173 CHOL & $\Sexpr{kf[3]}$ & $5.1 \times 10^7$ & $2.8 \times 10^{-4}$ & 38 & 0 & -1 & 1.21 \\
174 SM & $\Sexpr{kf[4]}$ & $3.7 \times 10^6$ & $3.1 \times 10^{-3}$ & 61 & 0 & 3 & 0.8 \\
175 PE & $\Sexpr{kf[5]}$ & $2.3 \times 10^6$ & $1 \times 10^{-5}$ & 55 & 0 & 0 & 1.33 \\
178 \caption{Kinetic parameters and molecular properties of lipid types}
179 \label{tab:kinetic_parameters_lipid_types}
182 %%% \DLA{I think we may just reduce these three sections; area, $k_\mathrm{f}$
183 %%% and $k_\mathrm{b}$ to \cref{tab:kinetic_parameters_lipid_types} with
186 %%% \RZ{Yes, but we also have to have then as comments the numbers that
187 %%% are not supported by refs}
189 \subsubsection{$k_\mathrm{f}$ for lipid types}
190 $k_{\mathrm{f}_\mathrm{PC}}$ was measured by
191 \citet{Nichols1985:phospholipid_monomer_vesicle_thermodynamics} and
192 was found to be $3.7\times 10^6$~$\frac{1}{\mathrm{M}\mathrm{s}}$ by the
193 partitioning of \ac{NBDPC} between \ac{DOPC} vesicles and water. As
194 similar references are not available for \ac{SM} or \ac{PS}, we assume that they have
195 the same $k_\mathrm{f}$. For \ac{CHOL}, no direct measurement of $k_\mathrm{f}$ is available,
196 however, \citet{Estronca2007:dhe_kinetics} measured the transfer of
197 \ac{DHE} from \ac{BSA} to \acp{LUV}, and found a $k_\mathrm{f}$ of $5.1\times 10^7$~%
198 $\frac{1}{\mathrm{M} \mathrm{s}}$. We assume that this value is close
199 to that of \ac{CHOL}, and use it for $k_{\mathrm{f}_\mathrm{\ac{CHOL}}}$. In the case of
200 \ac{PE}, \citet{Abreu2004:kinetics_ld_lo} measured the association of
201 \ac{NBDDMPE} with \ac{POPC} \acp{LUV} and found a value for $k_\mathrm{f}$ of $2.3 \times 10^{6}$~%
202 $\frac{1}{\mathrm{M} \mathrm{s}}$. These three authors used a slightly
203 different kinetic formalism ($\frac{d\left[A_v\right]}{dt} =
204 k'_\mathrm{f}[A_m][L_v] - k_\mathrm{b}[A_v]$), so we correct their values of $k_\mathrm{f}$ by
205 multiplying by the surface area of a mole of lipids.
207 \subsubsection{$k_\mathrm{b}$ for lipid types}
209 \citet{Wimley1990:dmpc_exchange} measured the half time for the
210 exchange of [$^3$H]\ac{DMPC} between \acp{LUV} at 30\textdegree C, and found it
211 to be 9.6 hr. As this is a first order reaction, and the primary
212 limiting step in exchange is lipid desorption, $k_\mathrm{b}$ for \ac{DMPC} is
213 $k_{\mathrm{b}_\mathrm{PC}}=\frac{\log 2}{9.6 \times 60 \times 60} \approx
214 \Sexpr{log(2)/(9.6*60*60)}
215 \mathrm{s}^{-1}$. We assume that $k_\mathrm{b}$ for \ac{SM} is the same as for \ac{PC}.
216 To estimate the $k_\mathrm{b}$ of \ac{PE} and \ac{PS}, we used the data from
217 \citet{Nichols1982:ret_amphiphile_transfer} who measured the rate of
218 exchange of the fluorescent label \ac{C6NBD} attached to different
219 lipid species. Although the values of $k_\mathrm{b}$ are different for the labeled
220 and unlabeled lipids, we assume that the ratios of the kinetics
221 constants for the lipid types are the same. Furthermore we assume that
222 \ac{PG} behaves similarly to \ac{PS}. Thus, we can determine the $k_\mathrm{b}$ of \ac{PE} and
223 \ac{PS} from the already known $k_\mathrm{b}$ of \ac{PC}. For \ac{C6NBD} labeled \ac{PC},
224 \citet{Nichols1982:ret_amphiphile_transfer} obtained a $k_\mathrm{b}$ of
225 $0.89$~$\mathrm{min}^{-1}$, \ac{PE} of $0.45$~$\mathrm{min}^{-1}$ and PG of
226 $0.55$~$\mathrm{min}^{-1}$. Assuming the ratios of the rate of
227 exchange are the same for unlabeled lipids and labeled lipids, we can
228 determine the $k_\mathrm{b}$ of \ac{PE} and \ac{PS} from the already known $k_\mathrm{b}$ of
229 \ac{PC}~\citep{Wimley1990:dmpc_exchange}. Calculating the ratio, this leads
230 us to $k_{\mathrm{b}_\mathrm{PE}} =
231 \frac{k_{\mathrm{b}_\mathrm{PC}}\times\mathrm{PE}}{\mathrm{PC}} \approx
232 \frac{2\times 10^{-5}\,\mathrm{s}^{-1} \times
233 0.45\,\mathrm{min}^{-1}}{0.89\,\mathrm{min}^{-1}} \approx
234 \Sexpr{log(2)/(9.6*60*60)*0.45/0.89}$~$\mathrm{s}^{-1}$
235 and likewise, $k_{\mathrm{b}_\mathrm{PS}}\approx
236 \Sexpr{log(2)/(9.6*60*60)*0.55/0.89}$~$\mathrm{s}^{-1}$.
237 The $k_\mathrm{b}$ of \ac{SM} was determined using the work of
238 \citet{Bai1997:lipid_movementbodipy}, who measured spontaneous
239 transfer of C$_5$-DMB-SM and C$_5$-DMB-PC from donor and acceptor
240 vesicles, finding $3.4\times10^{-2}$~$\mathrm{s}^{-1}$ and
241 $2.2\times10^{-3}$~$\mathrm{s}^{-1}$ respectively; using the ratio of
242 $k_\mathrm{b}$ of C$_5$-DMB-SM to the $k_\mathrm{b}$ of C$_5$-DMB-PC times the $k_\mathrm{b}$ of
243 \ac{PC} ($\frac{3.4 \times 10^{-2}\mathrm{s}^{-1}}{2.2 \times
244 10^{-3}\mathrm{s}^{-1}}\approx
245 \Sexpr{log(2)/(9.6*60*60)}\mathrm{s}^{-1}$),
246 we obtain $k_{\mathrm{b}_\mathrm{SM}} \approx
247 \Sexpr{log(2)/(9.6*60*60)*3.4e-2/2.2e-3}$.
248 In the case of \ac{CHOL}, \citet{Jones1990:lipid_transfer} measured the
249 $t_{1/2}$ of [$^3$H] transfer from POPC vesicles and found it to be 41
250 minutes, leading to a $k_{\mathrm{b}_\mathrm{CHOL}} = \frac{\log 2}{41\times
252 \Sexpr{log(2)/(41*60)}$~$\mathrm{s}^{-1}$.
254 \subsubsection{Headgroup Surface Area for lipid types}
256 % Smondyrev, A. M., and M. L. Berkowitz. 1999b. United atom force
257 % field for phospholipid membranes: constant pressure molecular
258 % dynamics simulation of dipalmitoylphosphatidicholine/water system.
259 % J. Comput. Chem. 50:531–545.
261 Different lipids have different headgroup surface areas, which contributes to
262 $\left[S_\mathrm{vesicle}\right]$. \citet{Smaby1997:pc_area_with_chol}
263 measured the surface area of \ac{POPC} with a Langmuir film balance, and
264 found it to be 63~Å$^2$ at $30$~$\frac{\mathrm{mN}}{\mathrm{m}}$.
265 Molecular dynamic simulations found an area of 54 Å$^2$ for
266 \ac{DPPS}\citep{Cascales1996:mds_dpps_area,Pandit2002:mds_dpps}, which is
267 in agreement with the experimental value of 56~Å$^2$ found using a
268 Langmuir balance by \citet{Demel1987:ps_area}.
269 \citet{Shaikh2002:pe_phase_sm_area} measured the area of \ac{SM} using a
270 Langmuir film balance, and found it to be 61~Å$^2$. Using $^2$H NMR,
271 \citet{Thurmond1991:area_of_pc_pe_2hnmr} found the area of
272 \ac{DPPE}-d$_{62}$ to be 55.4 Å$^2$. \citet{Robinson1995:mds_chol_area}
273 found an area for \ac{CHOL} of 38~Å$^2$ using molecular dynamic
276 % robinson's chol area is kind of crappy; they did it using MDS, but
277 % vaguely; they don't indicate what it was at the end, and they only
278 % had 2 molecules of \ac{CHOL}.
282 % \ac{PC}: 63 A (Smaby97rd) (Pandit (?)
283 % \ac{PS}: 54 A (Pandit02LIB) (Cascales 1996; J. Chem. Phys 104:2713, Mukhopadhyay2004)
284 % \ac{CHOL}: 38 A (Robinson 1995) (Lundberg 1982)
285 % \ac{SM}: 51 A (Shaikh2002, but 61(?))
286 % PE: 55 A (Thurmond, 1991) (Pandit2002)
288 % Compare to results by Petrache2004 with DOPC of 72.5, DOPS of 65.3.
290 % The cross‑sectional area of cholesterol is ~37 ֵ2 (Lundberg, B. 1982. A
291 % surface film study of the lateral packing of phosphatidylcholine and
292 % cholesterol. Chem. Phys. Lipids. 31:23-32).
294 % From Pandit02LIB: the area per headgroup for PE is ~10‑20\% (Thurmond
295 % et al., 1991) smaller than the area per PC molecule. For DLPE, it is
296 % 50.6 ֵ2 (McIntosh and Simon, 1986; Damodaran et al., 1992). The
297 % estimated area per molecule for DPPE is ~55.4 ֵ2 (Thurmond et al.,
300 % From Pandit02LIB: (from their simulations) average area per headgroup
301 % for DPPS molecules is 53.75 ± 0.10 ֵ2. The values inferred from the
302 % experiments are in the region between 45 and 55 ֵ2 (Cevc et al., 1981.
303 % Titration of the phase transition of phosphatidylserine bilayer
304 % membranes. Effect of pH, surface electrostatics, ion binding, and
305 % head‑group hydration. Biochemistry. 20:4955‑4965; | Demel et al.,
306 % 1987. Monolayer characteristics and thermal behavior of natural and
307 % synthetic phosphatidylserines. Biochemistry. 26:8659‑8665). Note that
308 % the average area per headgroup for DPPC bilayer obtained from
309 % simulations and measurements is ~62 ֵ A2. One would expect that the
310 % area per headgroup in case of DPPS molecules will be larger than the
311 % area per DPPC molecule, because the headgroups of DPPS are charged and
312 % therefore repel each other. Contrary to this expectation, the area per
313 % headgroup in DPPS is ~13\% smaller than that of DPPC. López Cascales
314 % et al. (1996. Molecular dynamics simulation of a charged biological
315 % membrane. J. Chem. Phys. 104:2713‑2720.) proposed that this reduction
316 % in the area per headgroup is due to a strong intermolecular
317 % coordination between DPPS molecules. MD gave area per POPS of 55 A2 at
318 % 300K (Mukhopadhyay04LIB).
320 % However, by 2H NMR and X-ray, Petrache04LIB determined the area of
321 % DOPS at 30 C to be 65.3 A2, considerably smaller than DOPC, determined
322 % by the same group in another paper to be 72.5 ֵ A2
323 % (Tristram‑Nagle00LIB).
325 % From Shaikh02: At 35 mN/m of surface pressure, the areas/molecule for
326 % POPE, SM, and CHOL were found to be 65.2, 61.1, and 38.8 Å2.
328 % In Lonchin99: POPC area 72 A2 ((a) Cornell, B. A.; Middlehurst, J.;
329 % Separovic, F. Biochim. Biophys. Acta1980, 598, 405−410), oleic acid
330 % area 32 A2 ((a) Small, D. M. In Handbook of Lipid Research; Hanahan,
331 % D. J., Ed.; Plenum Press: New York, 1986; Vol. 4, Chapter 9.
333 % (b) Jain, M. K. Introduction to Biological Membranes, 2nd ed.; John
334 % Wiley \& Sons: New York, 1988; Chapter 3). Cistola, D. P.; Atkinson,
335 % D.; Hamilton, J. A.; Small, D. M. Biochemistry1986, 25, 2804−2812.
337 % From Marsh96: The values obtained for the headgroup area per molecule,
338 % fitted as free parameters in Fig. 6, are APE = 58 A2 and APC = 78 A2,
339 % for DOPE and DOPC, respectively. These values are within the range
340 % obtained for the corresponding areas per molecule of
341 % phosphatidylethanolamines and phosphatidylcholines in the fluid
342 % lamellar phase (see, e.g., Marsh, 1990), and correspond to values for
343 % the lipid packing parameter of V/Al = 1.08 and 1.34 for DOPC and DOPE,
344 % respectively. The former value is consistent with DOPC alone forming
345 % hydrated lamellar structures rather than nonlamellar structures. The
346 % latter value lies at the lower end of the packing parameters obtained
347 % for phosphatidylethanolamines.
349 % From Kumar91: The head group area chosen was 71.7 A2 for PC and
350 % lysoPC, 42 A2 for PE, and 19 A2 for Chol (5, 14). The value selected
351 % for PC is in agreement with the 71 A2 value determined for PCs with
352 % saturated acyl chains (14). The hydrocarbon volume and the length
353 % chosen for Chol are 400 A3 and 17.25 A, respectively (5).
357 % From Sampaio05: Besides this work and our own earlier report on the
358 % association of NBD-DMPE with lipid bilayers (Abreu et al., 2004), we
359 % are aware of only one other report in the literature (Nichols, 1985)
360 % in which all the kinetic constants of lipid-derived amphiphile
361 % association with lipid bilayer membranes were experimentally measured.
362 % {indeed; everything is k- !!!; rz}
364 % From McLean84LIB: Although it is difficult to measure cmc values for
365 % the sparingly soluble lipids used in this study, estimates for
366 % lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds,
367 % 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution;
368 % Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6
369 % X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X
370 % 10-8 M for DMPC was estimated from the linear relationship between ln
371 % cmc and the number of carbons in the PC acyl chain by using data for n
372 % = 7, 8, 10, and 16 [summarized in Tanford (1980)].
374 % From Toyota08: Recently, several research groups have reported
375 % self-reproducing systems of giant vesicles that undergo a series of
376 % sequential transformations: autonomous growth, self-division, and
377 % chemical reactions to produce membrane constituents within the giant
380 % Vesicle sizes of 25 nm for SUV and 150 nm for LUV were mentioned by
383 % From Lund-Katz88: Charged and neutral small unilamellar vesicles
384 % composed of either saturated PC, unsaturated PC, or SM had similar
385 % size distributions with diameters of 23 \& 2 nm.
387 % From Sampaio05LIB: The exchange of lipids and lipid derivatives
388 % between lipid bilayer vesicles has been studied for at least the last
389 % 30 years. Most of this work has examined the exchange of amphiphilic
390 % molecules between a donor and an acceptor population. The measured
391 % efflux rates were shown in almost all cases, not surprisingly, to be
392 % first order processes. In all of this work, the insertion rate has
393 % been assumed to be much faster than the efflux rate. Having measured
394 % both the insertion and desorption rate constants for amphiphile
395 % association with membranes, our results show that this assumption is
396 % valid. In several cases reported in the literature, the insertion rate
397 % constant was assumed, although never demonstrated, to be a
398 % diffusion-controlled process.
400 % (for methods? From McLean84LIB: The activation free energies and free
401 % energies of transfer from self-micelles to water increase by 2.2 and
402 % 2.1 kJ mol-' per methylene group, respectively. {see if we can use it
403 % to justify arranging our changed in activating energy around 1
406 % Jones90 give diameter of LUV as 100 nm, and of SUV as 20 nm; that
407 % would give the number of molecules per outer leaflet of a vesicle as
410 % Form Simard08: Parallel studies with SUV and LUV revealed that
411 % although membrane curvature does have a small effect on the absolute
412 % rates of FA transfer between vesicles, the ΔG of membrane desorption
413 % is unchanged, suggesting that the physical chemical properties which
414 % govern FA desorption are dependent on the dissociating molecule rather
415 % than on membrane curvature. However, disagreements on this fundamental
416 % issue continue (57, 61, 63, 64)
418 % (methods regarding the curvature effect: Kleinfeld93 showed that the
419 % transfer parameters of long-chain FFA between the lipid vesicles
420 % depend on vesicle curvature and composition. Transfer of stearic acid
421 % is much slower from LUV as compared to SUV).
423 % From McLean84: In a well-defined experimental system consisting of
424 % unilamellar lipid vesicles, in the absence of protein, the
425 % rate-limiting step for the overall exchange process is desorption
426 % (McLean \& Phillips, 1981). {thus I can take exchange data for the
427 % estimation of k- rz; 8/11/08}.
429 \subsubsection{Complex Formation 1}
431 \citet{Thomas1988:chol_transfer} found that \ac{SM} significantly decreases
432 the rate of cholesterol transfer from \ac{PC}, while \ac{PE} and \ac{PS} have no
433 effect at physiologically significant levels. We attribute this to the
434 formation of complexes between \ac{SM} and \ac{CHOL}, which retards the rate of
435 efflux of \ac{CHOL} from the membrane. Similar complexes exist between \ac{PC}
436 and \ac{CHOL}, where it was shown that two \ac{CHOL} molecules complex with one
437 \ac{PC}~\citep{Huang1999:chol_solubility_model,
438 Huang1999:cholesterol_solubility,McConnell2006,McConnell2003}. It
439 has also been shown that \ac{SM} binds more \ac{CHOL} molecules than
440 \ac{PC}~\citep{Katz1988:pl_packing_chol}. We assume that three \ac{CHOL} complex
441 with one \ac{SM}, leading to $\mathrm{CF}1$ values of $3$, $2$, and $-1$
442 for \ac{SM}, \ac{PC}, and \ac{CHOL}, respectively.
444 \subsubsection{Curvature}
446 We used the formulation for curvature of
447 \citet{Israelachvili1975:amphiphile_self_assembly}, or $S=\frac{v}{l_c
448 a}$, where $S$ is the curvature, which ranges from $0$ to $\infty$,
449 $l_c$ is the critical length of the acyl chain, $v$ is the volume of
450 the lipid, and $a$ is the area of the head group.
451 \citet{Kumar1991:lipid_packing} found a value for $S$ of 1.33 for PE,
452 1.21 for \ac{CHOL}, and 0.8 for diC$_{16}$ \ac{PC}. We assume that \ac{PS} has neutral
453 curvature (value of 1), and that \ac{SM} has the same curvature as \ac{PC}
454 (0.8). It should also be noted that in reality the curvature parameter
455 changes with length, but at longer chain lengths, is effectively
456 constant; in the current model, curvature is held constant for each
460 % 1. Israelachvili, J. N., Mitchell, D. J. & Ninham, B. W. (1976) J.
461 % Chem. Soc. Faraday Trans. 2 72, 1525-1568.
462 % 2. Israelachvili, J. N., Marcelja, S. & Horn, R. G. (1980) Q. Rev.
463 % Biophys. 13, 121-200.
464 % 3. Israelachvili, J. N. (1985) Intermolecular and Surface Forces
465 % (Academic, New York), pp. 229-271.
467 % \DLA{This section needs to be added still.}
469 % PE = 1.33 (Kumar91)
470 % CHOL = 1.21 (Kumar91)
472 % SM=0.8 (assumed by rz same as PC)
473 % PS=1 (no refs so far; should be close to unity; rz)
475 % From Epand94LIB: a considerable amount of heat is evolved per mole of
476 % bilayer stabilizer when such molecules are inserted into membranes
477 % which are under a high curvature strain. If this energy were coupled
478 % to a membrane event, such as the conformational change in a protein,
479 % it could be the driving force responsible for such processes.
481 % From Epand94LIB: considerable energy may be released upon the
482 % incorporation of certain molecules into membranes which have a low
483 % radius of spontaneous curvature. {discuss with DL in terms of
484 % catalytic activity; rz; 8/24/2010}.
486 % Kleinfeld93 showed that the transfer parameters of long-chain FFA
487 % between the lipid vesicles depend on vesicle curvature and
488 % composition. Transfer of stearic acid is much slower from LUV as
491 % Substitution of DPPC SUV for POPC SUV as the donor matrix revealed
492 % some differences in the energetics of transfer that are associated
493 % with the physical state of the donor particle. Whereas the free
494 % energies of transfer of oleic acid from DPPC SUV and POPC SUV are
495 % similar, the free energy of transfer of stearic acid from DPPC SUV is
496 % much higher than that observed from POPC.
498 % From Kumar91: S is given by S = V/al, where V is the hydrocarbon
499 % volume, a is the area of head group, and l is the critical length of
500 % the hydrocarbon chain (1-3). a, V, and I are all estimable or
501 % measurable (4), and the value of S can be calculated. The value of S
502 % determines the aggregate formed by lipids or any amphiphiles upon
503 % hydration. It has been shown that lipids aggregate to form spherical
504 % micelles (S < 1/3), nonspherical (cylindrical) micelles (1/3 < S <
505 % 1/2), bilayers (1/2 < S < 1), and reverse micelles or hexagonal (HII,)
506 % phases (S > 1). However, theoreticians caution that the above
507 % predicted limits set on the values of S are relatively insensitive to
508 % the exact values of V and a but are strongly dependent upon the choice
511 % From Kumar91: Cylindrical and wedge-shaped molecules have been shown
512 % to be essential for spontaneous vesiculation (6) as well as bilayer
513 % stabilization (7, 8).
515 % From Kumar91: It has also been shown experimentally that molecules
516 % having complementary molecular shapes can form bilayer structures.
517 % Inverted cone (wedge)- shaped lysoPC and cone-shaped Chol or
518 % unsaturated PE can interact stoichiometrically to form bilayer
519 % structures (9-12). {quote to justify our formula for complementary
522 \section{Kinetic adjustments}
523 \label{sec:kinetic_adjustments}
525 In the mass action equation used in the formalism, both the forward
526 and backwards kinetic constants ($k_\mathrm{f}$ and $k_\mathrm{b}$) are adjusted (by
527 $k_{\mathrm{f}i\mathrm{adj}}$ and $k_{\mathrm{b}i\mathrm{adj}}$) to reflect the
528 properties of the vesicle and the properties of the monomer entering
529 or exiting the vesicle.
531 \subsection{Forward Kinetic Adjustments}
533 The forward rate constant adjustment, $k_{\mathrm{f}i\mathrm{adj}}$ takes into
534 account unsaturation ($un_\mathrm{f}$), charge ($ch_\mathrm{f}$), curvature ($cu_\mathrm{f}$),
535 length ($l_\mathrm{f}$), and complex formation ($CF1_\mathrm{f}$), each of which is
536 modified depending on the specific component and the vesicle into
537 which the component is entering; the final $k_{\mathrm{f}i\mathrm{adj}}$ is the
538 product of the adjustments made separately for each property.
541 k_{\mathrm{f}i\mathrm{adj}} = un_\mathrm{f} \cdot ch_\mathrm{f} \cdot cu_\mathrm{f} \cdot l_\mathrm{f} \cdot CF1_\mathrm{f}
545 \subsubsection{Unsaturation Forward}
547 In order for a lipid to be inserted into a membrane, a void has to be
548 formed for it to fill. Voids can be generated by the combination of
549 unsaturated and saturated lipids forming heterogeneous domains. Void
550 formation is increased when the unsaturation of lipids in the vesicle
551 is widely distributed; in other words, the insertion of lipids into
552 the membrane is greater when the standard deviation of the
553 unsaturation is larger. %
554 %%% \RZ{May need to site (at least for us) works showing
555 %%% mismatch-dependent ``defects''}. %
556 Assuming that an increase in width of the distribution linearly
557 decreases the free energy of activation, the $un_\mathrm{f}$ parameter must
558 follow $a^{\mathrm{stdev}\left(un_\mathrm{vesicle}\right)}$ where $a > 1$,
559 so a convenient starting base for $a$ which leads to a reasonable
560 change in Eyring activation energy is $2$:
563 un_\mathrm{f} = 2^{\mathrm{stdev}\left(un_\mathrm{vesicle}\right)}
564 \label{eq:unsaturation_forward}
567 The mean $\mathrm{stdev}\left(un_\mathrm{vesicle}\right)$ in our
568 simulations is around $1.5$, which leads to a $\Delta \Delta
569 G^\ddagger$ of $\Sexpr{to.latex(format(digits=3,to.kcal(2^1.5)))}
570 \frac{\mathrm{kcal}}{\mathrm{mol}}$, and a total range of possible
571 values depicted in \cref{fig:unf_graph}.
573 % \RZ{Explain here, or even earlier that the formulas were ad hoc
574 % adjusted to correspond to ``reasonable'' changes in the Eyring
575 % activation energies.}
577 % It is not clear that the unsaturation of the inserted monomer will
578 % affect the rate of the insertion positively or negatively, so we do
579 % not include a term for it in this formalism.
583 <<unf_graph,echo=FALSE,results='hide',fig.width=10,fig.height=5,cache=FALSE>>=
584 layout(matrix(1:2,nrow=1,ncol=2))
585 curve(2^x,from=0,to=sd(c(0,4)),
586 xlab=expression(paste(stdev,group("(",italic(un[vesicle]),")"))),
587 ylab=expression(italic(un[f])))
588 vps <- baseViewports()
589 pushViewport(vps$figure)
591 just=c("left","top"),
592 gp=gpar(fontsize=48),
593 y=unit(1,"npc")-unit(0.25,"lines"),
594 x=unit(0,"npc")+unit(0.25,"lines"))
596 curve(to.kcal(2^x),from=0,to=sd(c(0,4)),
597 xlab=expression(paste(stdev,group("(",italic(un[vesicle]),")"))),
598 ylab="Activation Energy (kcal/mol)")
599 vps <- baseViewports()
600 pushViewport(vps$figure)
602 just=c("left","top"),
603 gp=gpar(fontsize=48),
604 y=unit(1,"npc")-unit(0.25,"lines"),
605 x=unit(0,"npc")+unit(0.25,"lines"))
608 \caption{Change in unsaturation forward ($un_\mathrm{f}$) (A) or activation
609 energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the
610 standard deviation of the unsaturation of the vesicle
611 ($\mathrm{stdev}\left(un_\mathrm{vesicle}\right)$).}
612 \label{fig:unf_graph}
615 \subsubsection{Charge Forward}
617 A charged lipid such as \ac{PS} approaching a vesicle with an average
618 charge of the same sign will experience repulsion, whereas those with
619 different signs will experience attraction, the degree of which is
620 dependent upon the charge of the monomer and the average charge of the
621 vesicle. If either the vesicle or the monomer has no charge, there
622 should be no effect of charge upon the rate. This leads us to the
623 following equation, $a^{-\left<ch_v\right> ch_m}$ (which is similar to
624 the numerator of Coulomb's law), where $\left<ch_v\right>$ is the
625 average charge of the vesicle, and $ch_m$ is the charge of the
626 monomer. If either $\left<ch_v\right>$ or $ch_m$ is 0, the adjustment
627 parameter is 1 (no change), whereas it decreases if both are positive
628 or negative, as the product of two real numbers with the same sign is
629 always positive. The base for $a$ ($60$) was chosen ad hoc to
630 correspond to a maximum of about
631 $0.5\frac{\mathrm{kcal}}{\mathrm{mol}}$ change in activation energy
632 for the common case, resulting in the following equation:
636 ch_\mathrm{f} = 60^{-\left<{ch}_v\right> {ch}_m}
637 \label{eq:charge_forward}
640 The most common $\left<{ch}_v\right>$ is around $-0.165$, which leads
641 to a range of $\Delta \Delta G^\ddagger$ from
642 $\Sexpr{format(digits=3,to.kcal(60^(-.165*-1)))}
643 \frac{\mathrm{kcal}}{\mathrm{mol}}$ to
644 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$, and the total range of possible
645 values seen in \cref{fig:chf_graph}.
649 <<chf_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
650 trellis.device(new=F,color=TRUE)
651 trellis.par.set(list(axis.line =list(col="transparent")))
652 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
653 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
654 x <- seq(-1,0,length.out=20)
655 y <- seq(-1,0,length.out=20)
656 grid <- expand.grid(x=x,y=y)
657 grid$z <- as.vector(60^(-outer(x,y)))
658 print(wireframe(z~x*y,grid,cuts=50,
661 scales=list(arrows=F,col=1,distance=0.75),
662 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
663 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
664 zlab=list(label=expression(italic(ch[f])),rot=93)),newpage=FALSE)
667 just=c("left","top"),
668 gp=gpar(fontsize=48),
669 y=unit(1,"npc")-unit(0.25,"lines"),
670 x=unit(0,"npc")+unit(0.25,"lines"))
672 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
673 x <- seq(-1,0,length.out=20)
674 y <- seq(-1,0,length.out=20)
675 grid <- expand.grid(x=x,y=y)
676 grid$z <- as.vector(to.kcal(60^(-outer(x,y))))
677 print(wireframe(z~x*y,grid,cuts=50,
680 scales=list(arrows=F,col=1,distance=0.75),
681 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
682 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
683 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
686 just=c("left","top"),
687 gp=gpar(fontsize=48),
688 y=unit(1,"npc")-unit(0.25,"lines"),
689 x=unit(0,"npc")+unit(0.25,"lines"))
692 \caption{Change in charge forward ($ch_\mathrm{f}$) (A) or activation energy in
693 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the charge of the
694 monomer entering ($ch_\mathrm{monomer}$) and the average charge of
695 the vesicle ($\left<ch_\mathrm{vesicle}\right>$).}
696 \label{fig:chf_graph}
699 \subsubsection{Curvature Forward}
701 Curvature is a measure of the intrinsic propensity of specific lipids
702 to form micelles (positive curvature), inverted micelles (negative
703 curvature), or planar sheets (neutral
704 curvature)~\citep{Israelachvili1975:amphiphile_self_assembly}. %
705 In this formalism, curvature is measured as the ratio of the volume of
706 the lipid to the area of the head times the length of the lipid
707 ($S=\frac{v}{l_c a}$), so negative curvature is bounded by $(0,1)$,
708 neutral curvature is 1, and positive curvature is bounded by
709 $(1,\infty)$. The curvature can be transformed using $\log$, which has
710 the property of making the range of positive and negative curvature
711 equal, and distributed about 0.
713 As in the case of unsaturation, void formation is increased by the
714 presence of lipids with mismatched curvature. Thus, a larger
715 distribution of curvature in the vesicle increases the rate of lipid
716 insertion into the vesicle. However, a component with curvature
717 $cu^{-1}$ will cancel out a component with curvature $cu$, so we have
718 to log transform (turning these into $-\log cu$ and $\log cu$), then
719 take the absolute value ($\log cu$ and $\log cu$), and finally measure
720 the width of the distribution, which in the case of exactly mismatched
721 curvatures is $0$. Thus, by using the log transform to make the range
722 of the lipid curvature equal between positive and negative, and taking
723 the average to cancel out exactly mismatched curvatures, we come to an
724 equation with the shape $a^{\left<\log cu_\mathrm{vesicle}\right>}$.
725 An ad hoc base for $a$ which for the common case is $10$, yielding:
729 % cu_\mathrm{f} = 10^{\mathrm{stdev}\left|\log cu_\mathrm{vesicle}\right|}
730 cu_\mathrm{f} = 10^{\left|\left<\log cu_\mathrm{monomer} \right>\right|\mathrm{stdev} \left|\log cu_\mathrm{vesicle}\right|}
731 \label{eq:curvature_forward}
734 The most common $\left|\left<\log {cu}_v\right>\right|$ is around
735 $0.013$, which with the most common $\mathrm{stdev} \log
736 cu_\mathrm{vesicle}$ of $0.213$ leads to a $\Delta \Delta G^\ddagger$
737 of $\Sexpr{format(digits=3,to.kcal(10^(0.13*0.213)))}
738 \frac{\mathrm{kcal}}{\mathrm{mol}}$. This is a consequence of the
739 relatively matched curvatures in our environment. The full range of
740 $cu_\mathrm{f}$ values possible are shown in \cref{fig:cuf_graph}.
742 % 1.5 to 0.75 3 to 0.33
744 <<cuf_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
745 trellis.device(new=F,color=TRUE)
746 trellis.par.set(list(axis.line =list(col="transparent")))
747 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
748 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
749 grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
750 sd(abs(log(c(1,0.33)))),
751 sd(abs(log(c(0.33,3)))))),length.out=20),
752 y=seq(0,max(c(mean(log(c(1,3)),
753 mean(log(c(1,0.33))),
754 mean(log(c(0.33,3)))))),length.out=20))
755 grid$z <- 10^(grid$x*grid$y)
756 print(wireframe(z~x*y,grid,cuts=50,
759 scales=list(arrows=F,col=1,distance=0.75),
760 xlab=list(label=expression(paste(stdev,group("(",italic(log(cu[vesicle])),")"))),rot=30),
761 ylab=list(label=expression(paste(symbol("\341"),italic(log(cu[vesicle])),symbol("\361"))),rot=-35),
762 zlab=list(label=expression(italic(cu[f])),rot=93)),newpage=FALSE)
765 just=c("left","top"),
766 gp=gpar(fontsize=48),
767 y=unit(1,"npc")-unit(0.25,"lines"),
768 x=unit(0,"npc")+unit(0.25,"lines"))
770 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
771 grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
772 sd(abs(log(c(1,0.33)))),sd(abs(log(c(0.33,3)))))),length.out=20),
773 y=seq(0,max(c(mean(log(c(1,3)),
774 mean(log(c(1,0.33))),
775 mean(log(c(0.33,3)))))),length.out=20))
776 grid$z <- to.kcal(10^(grid$x*grid$y))
777 print(wireframe(z~x*y,grid,cuts=50,
780 scales=list(arrows=F,col=1,distance=0.75),
781 xlab=list(label=expression(paste(stdev,group("(",italic(log(cu[vesicle])),")"))),rot=30),
782 ylab=list(label=expression(paste(symbol("\341"),italic(log(cu[vesicle])),symbol("\361"))),rot=-35),
783 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
786 just=c("left","top"),
787 gp=gpar(fontsize=48),
788 y=unit(1,"npc")-unit(0.25,"lines"),
789 x=unit(0,"npc")+unit(0.25,"lines"))
792 \caption{Change in curvature forward ($cu_\mathrm{f}$) (A) or activation energy
793 in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the standard
794 deviation of the log of the curvature of the vesicle
795 ($\mathrm{stdev}\left(\log cu_\mathrm{vesicle}\right)$) and the mean
796 of the log of the curvature of the vesicle ($\left<\log
797 cu_\mathrm{vesicle}\right>$).}
798 \label{fig:cuf_graph}
801 \subsubsection{Length Forward}
803 As in the case of unsaturation, void formation is easier when vesicles
804 are made up of components of mismatched lengths. Thus, when the width
805 of the distribution of lengths is larger, the forward rate should be
806 greater as well, leading us to an equation of the form
807 $x^{\mathrm{stdev} l_\mathrm{vesicle}}$, where $\mathrm{stdev}
808 l_\mathrm{vesicle}$ is the standard deviation of the length of the
809 components of the vesicle, which has a maximum possible value of
810 $\Sexpr{format(digits=3,sd(c(rep(12,50),rep(24,50))))}$ and a minimum
811 of 0 in this set of simulations. Based on activation energy, a
812 reasonable base for $x$ is 2, leading to:
815 l_\mathrm{f} = 2^{\mathrm{stdev} l_\mathrm{vesicle}}
816 \label{eq:length_forward}
819 The most common $\mathrm{stdev} l_\mathrm{vesicle}$ is around $3.4$, which leads to
820 a $\Delta \Delta G^\ddagger$ of
821 $\Sexpr{format(digits=3,to.kcal(2^(3.4)))}
822 \frac{\mathrm{kcal}}{\mathrm{mol}}$.
824 % While it could be argued that increased length of the monomer could
825 % affect the rate of insertion into the membrane, it's not clear whether
826 % it would increase (by decreasing the number of available hydrogen
827 % bonds, for example) or decrease (by increasing the time taken to fully
828 % insert the acyl chain, for example) the rate of insertion or to what
829 % degree, so we do not take it into account in this formalism.
831 % \DLA{Incorporate McLean84 here}
832 % From McLean84LIB: Although it is difficult to measure cmc values for
833 % the sparingly soluble lipids used in this study, estimates for
834 % lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds,
835 % 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution;
836 % Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6
837 % X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X
838 % 10-8 M for DMPC was estimated from the linear relationship between ln
839 % cmc and the number of carbons in the PC acyl chain by using data for n
840 % = 7, 8, 10, and 16 [summarized in Tanford (1980)].
842 % From Nichols85: The magnitude of the dissociation rate constant
843 % decreases by a factor of approximately 3.2 per carbon increase in acyl
844 % chain length (see Table II here) {take into account for the formula;
847 % From Nichols85: The magnitude of the partition coefficient increases
848 % with acyl chain length [Keq(M-C6-NBD-PC) = (9.8±2.1) X 106 M-1 and Keq
849 % (P-C6-NBD-PC) = (9.4±1.3) x 107 M-1. {take into account for the
850 % formula; rz 8/17/2010}.
852 % From Nichols85: The association rate constant is independent of acyl
853 % chain length. {take into account for the formula; rz 8/17/2010}.
858 <<lf_graph,echo=FALSE,results='hide',fig.width=14,fig.height=5>>=
859 layout(matrix(1:2,nrow=1,ncol=2))
860 curve(2^x,from=0,to=sd(c(12,24)),
861 xlab=expression(paste(stdev,group("(",italic(l[vesicle]),")"))),
862 ylab=expression(italic(l[f])),
865 vps <- baseViewports()
866 pushViewport(vps$figure)
868 just=c("left","top"),
869 gp=gpar(fontsize=48),
870 y=unit(1,"npc")-unit(0.25,"lines"),
871 x=unit(0,"npc")+unit(0.25,"lines"))
873 curve(to.kcal(2^x),from=0,to=sd(c(12,24)),
874 xlab=expression(paste(stdev,group("(",italic(l[vesicle]),")"))),
875 ylab="Activation Energy (kcal/mol)",
877 vps <- baseViewports()
878 pushViewport(vps$figure)
880 just=c("left","top"),
881 gp=gpar(fontsize=48),
882 y=unit(1,"npc")-unit(0.25,"lines"),
883 x=unit(0,"npc")+unit(0.25,"lines"))
886 \caption{Change in length forward ($l_\mathrm{f}$) (A) or activation energy in
887 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the standard
888 deviation of the length of the vesicle
889 ($\mathrm{stdev}\left(l_\mathrm{vesicle}\right)$).}
893 \subsubsection{Complex Formation}
894 There is no contribution of complex formation to the forward reaction
895 rate in the current formalism.
899 \label{eq:complex_formation_forward}
902 \subsection{Backward adjustments ($k_{\mathrm{b}i\mathrm{adj}}$)}
904 Just as the forward rate constant adjustment $k_{\mathrm{f}i\mathrm{adj}}$
905 does, the backwards rate constant adjustment $k_{\mathrm{b}i\mathrm{adj}}$
906 takes into account unsaturation ($un_\mathrm{b}$), charge ($ch_\mathrm{b}$), curvature
907 ($cu_\mathrm{b}$), length ($l_\mathrm{b}$), and complex formation ($CF1_\mathrm{b}$), each of
908 which is modified depending on the specific component and the vesicle
909 from which the component is exiting:
913 k_{\mathrm{b}i\mathrm{adj}} = un_\mathrm{b} \cdot ch_\mathrm{b} \cdot cu_\mathrm{b} \cdot l_\mathrm{b} \cdot CF1_\mathrm{b}
917 \subsubsection{Unsaturation Backward}
919 Unsaturation also influences the ability of a lipid molecule to leave
920 a membrane. If a molecule has an unsaturation level which is different
921 from the surrounding membrane, it will be more likely to leave the
922 membrane. The more different the unsaturation level is, the greater
923 the propensity for the lipid molecule to leave. However, a vesicle
924 with some unsaturation (eg. average unsaturation of 2) is more
925 favorable for lipids with more unsaturation (eg. 3) than lipids with
926 equivalently less unsaturation (eg. 1), so the difference in energy
927 between unsaturation is not linear. Therefore, an equation with the
928 shape $x^{\left| y^{-\left< un_\mathrm{vesicle}\right>
929 }-y^{-un_\mathrm{monomer}} \right| }$, where
930 $\left<un_\mathrm{vesicle}\right>$ is the average unsaturation of the
931 vesicle and $un_\mathrm{monomer}$ is the unsaturation of the monomer,
932 allows for increasing the efflux of molecules from membranes where
933 they strongly mismatch, while allowing vesicles with greater
934 unsaturation to tolerate greater unsaturation mismatch between
935 monomers and the vesicle. The bases x and y were chosen ad hoc to
936 produce reasonable Eyring energies of activation over the range of
937 unsaturations expected, leading to:
941 un_\mathrm{b} = 7^{1-\left(20\left(2^{-\left<un_\mathrm{vesicle} \right>} - 2^{-un_\mathrm{monomer}} \right)^2+1\right)^{-1}}
942 \label{eq:unsaturation_backward}
945 The most common $\left<un_\mathrm{vesicle}\right>$ is around $1.7$, which
946 leads to a range of $\Delta \Delta G^\ddagger$ from
947 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-0)^2+1))))}
948 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with 0 unsaturation
950 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-4)^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
951 for monomers with 4 unsaturations. See \cref{fig:unb_graph} for
952 the full range of possible values.
956 <<unb_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
957 trellis.device(new=F,color=TRUE)
958 trellis.par.set(list(axis.line =list(col="transparent")))
959 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
960 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
961 grid <- expand.grid(x=seq(0,4,length.out=20),
962 y=seq(0,4,length.out=20))
963 grid$z <- (7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1)))
964 print(wireframe(z~x*y,grid,cuts=50,
967 scales=list(arrows=F,col=1,distance=0.75),
968 xlab=list(label=expression(paste(symbol("\341"),italic(un[vesicle]),symbol("\361"))),rot=30),
969 ylab=list(label=expression(italic(un[monomer])),rot=-35),
970 zlab=list(label=expression(italic(un[b])),rot=93)),newpage=FALSE)
972 just=c("left","top"),
973 gp=gpar(fontsize=48),
974 y=unit(1,"npc")-unit(0.25,"lines"),
975 x=unit(0,"npc")+unit(0.25,"lines"))
977 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
978 grid <- expand.grid(x=seq(0,4,length.out=20),
979 y=seq(0,4,length.out=20))
980 grid$z <- to.kcal((7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1))))
981 print(wireframe(z~x*y,grid,cuts=50,
984 scales=list(arrows=F,col=1,distance=0.75),
985 xlab=list(label=expression(paste(symbol("\341"),italic(un[vesicle]),symbol("\361"))),rot=30),
986 ylab=list(label=expression(italic(un[monomer])),rot=-35),
987 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
990 just=c("left","top"),
991 gp=gpar(fontsize=48),
992 y=unit(1,"npc")-unit(0.25,"lines"),
993 x=unit(0,"npc")+unit(0.25,"lines"))
996 \caption{Change in unsaturation backward ($un_\mathrm{b}$) (A) or activation energy in
997 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the unsaturation of the
998 monomer leaving ($un_\mathrm{monomer}$) and the average unsaturation of
999 the vesicle ($\left<un_\mathrm{vesicle}\right>$).}
1000 \label{fig:unb_graph}
1003 \subsubsection{Charge Backwards}
1004 As in the case of monomers entering a vesicle, opposites attract.
1005 Monomers leaving a vesicle leave faster if their charge has the same
1006 sign as the average charge vesicle. An equation of the form
1007 $ch_\mathrm{b} = a^{\left<ch_v\right> ch_m}$ is then appropriate, and
1008 using a base of $a=20$ yields:
1011 ch_\mathrm{b} = 20^{\left<{ch}_v\right> {ch}_m}
1012 \label{eq:charge_backwards}
1015 The most common $\left<ch_v\right>$ is around $-0.164$, which leads to
1016 a range of $\Delta \Delta G^\ddagger$ from
1017 $\Sexpr{format(digits=3,to.kcal(20^(-.164*-1)))}
1018 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with charge $-1$ to
1019 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with charge $0$.
1020 See \cref{fig:chb_graph} for the full range of possible values
1025 <<chb_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1026 trellis.device(new=F,color=TRUE)
1027 trellis.par.set(list(axis.line =list(col="transparent")))
1028 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1029 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1030 x <- seq(-1,0,length.out=20)
1031 y <- seq(-1,0,length.out=20)
1032 grid <- expand.grid(x=x,y=y)
1033 grid$z <- as.vector(20^(outer(x,y)))
1034 print(wireframe(z~x*y,grid,cuts=50,
1037 scales=list(arrows=F,col=1,distance=0.75),
1038 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
1039 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
1040 zlab=list(label=expression(italic(ch[b])),rot=93)),newpage=FALSE)
1043 just=c("left","top"),
1044 gp=gpar(fontsize=48),
1045 y=unit(1,"npc")-unit(0.25,"lines"),
1046 x=unit(0,"npc")+unit(0.25,"lines"))
1048 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1049 x <- seq(-1,0,length.out=20)
1050 y <- seq(-1,0,length.out=20)
1051 grid <- expand.grid(x=x,y=y)
1052 grid$z <- to.kcal(as.vector(20^(outer(x,y))))
1053 print(wireframe(z~x*y,grid,cuts=50,
1056 scales=list(arrows=F,col=1,distance=0.75),
1057 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
1058 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
1059 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1062 just=c("left","top"),
1063 gp=gpar(fontsize=48),
1064 y=unit(1,"npc")-unit(0.25,"lines"),
1065 x=unit(0,"npc")+unit(0.25,"lines"))
1068 \caption{Change in charge backward ($ch_\mathrm{b}$) (A) or activation energy in
1069 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the charge of the
1070 monomer leaving ($ch_\mathrm{monomer}$) and the average charge of
1071 the vesicle ($\left<ch_\mathrm{vesicle}\right>$).}
1072 \label{fig:chb_graph}
1076 \subsubsection{Curvature Backwards}
1078 The less a monomer's intrinsic curvature matches the average curvature
1079 of the vesicle in which it is in, the greater its rate of efflux. If
1080 the curvatures match exactly, $cu_\mathrm{f}$ needs to be one. To map negative and
1081 positive curvature to the same range, we also need take the logarithm.
1082 Positive ($cu > 1$) and negative ($0 < cu < 1$) curvature lipids cancel each other out,
1083 resulting in an average curvature of 1; the average of the log also
1084 has this property (average curvature of 0). Increasing mismatches in curvature increase the
1085 rate of efflux, but asymptotically. An equation which satisfies these
1086 criteria has the form $cu_\mathrm{f} = a^{1-\left(b\left( \left< \log
1087 cu_\mathrm{vesicle} \right> -\log
1088 cu_\mathrm{monomer}\right)^2+1\right)^{-1}}$. An alternative
1089 form would use the absolute value of the difference, however, this
1090 yields a cusp and sharp increase about the curvature equilibrium. We
1091 have chosen bases of $a=7$ and $b=20$ ad hoc, based on activation
1092 energy considerations.
1095 cu_\mathrm{b} = 7^{1-\left(20\left(\left< \log cu_\mathrm{vesicle} \right> -\log cu_\mathrm{monomer} \right)^2+1\right)^{-1}}
1096 \label{eq:curvature_backwards}
1099 The most common $\left<\log cu_\mathrm{vesicle}\right>$ is around
1100 $-0.013$, which leads to a range of $\Delta \Delta G^\ddagger$ from
1101 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(0.8))^2+1))))}
1102 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature 0.8 to
1103 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature
1105 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(1.3))^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
1106 for monomers with curvature 1.3. The full range of values possible for
1107 $cu_\mathrm{b}$ are shown in \cref{fig:cub_graph}.
1109 % \RZ{What about the opposite curvatures that actually do fit to each
1113 % From Kumar91: It has also been shown experimentally that molecules
1114 % having complementary molecular shapes can form bilayer structures.
1115 % Inverted cone (wedge)- shaped lysoPC and cone-shaped Chol or
1116 % unsaturated PE can interact stoichiometrically to form bilayer
1117 % structures (9-12). {quote to justify our formula for complementary
1122 <<cub_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1123 trellis.device(new=F,color=TRUE)
1124 trellis.par.set(list(axis.line =list(col="transparent")))
1125 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1126 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1127 grid <- expand.grid(x=log(seq(0.8,1.33,length.out=20)),
1128 y=seq(0.8,1.33,length.out=20))
1129 grid$z <- 7^(1-1/(20*(grid$x-log(grid$y))^2+1))
1130 print(wireframe(z~x*y,grid,cuts=50,
1133 scales=list(arrows=F,col=1,distance=0.75),
1134 xlab=list(label=expression(paste(symbol("\341"),log(italic(cu[vesicle])),symbol("\361"))),rot=30),
1135 ylab=list(label=expression(italic(cu[monomer])),rot=-35),
1136 zlab=list(label=expression(italic(cu[b])),rot=93)),newpage=FALSE)
1139 just=c("left","top"),
1140 gp=gpar(fontsize=48),
1141 y=unit(1,"npc")-unit(0.25,"lines"),
1142 x=unit(0,"npc")+unit(0.25,"lines"))
1144 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1145 grid <- expand.grid(x=log(seq(0.8,1.33,length.out=20)),
1146 y=seq(0.8,1.33,length.out=20))
1147 grid$z <- to.kcal(7^(1-1/(20*(grid$x-log(grid$y))^2+1)))
1148 print(wireframe(z~x*y,grid,cuts=50,
1151 scales=list(arrows=F,col=1,distance=0.75),
1152 xlab=list(label=expression(paste(symbol("\341"),log(italic(cu[vesicle])),symbol("\361"))),rot=30),
1153 ylab=list(label=expression(italic(cu[monomer])),rot=-35),
1154 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1157 just=c("left","top"),
1158 gp=gpar(fontsize=48),
1159 y=unit(1,"npc")-unit(0.25,"lines"),
1160 x=unit(0,"npc")+unit(0.25,"lines"))
1163 \caption{Change in curvature backward ($cu_\mathrm{b}$) (A) or activation
1164 energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the
1165 curvature of the monomer leaving ($ch_\mathrm{monomer}$) and the
1166 average of the log of the curvature of the vesicle
1167 ($\left<\log cu_\mathrm{vesicle}\right>$).}
1168 \label{fig:cub_graph}
1172 \subsubsection{Length Backwards}
1174 In a model membrane, the dissociation constant increases by a factor
1175 of approximately 3.2 per carbon decrease in acyl chain length~%
1176 \citep{Nichols1985:phospholipid_monomer_vesicle_thermodynamics}.
1177 Unfortunately, the experimental data known to us only measures the
1178 effect of chain lengths less than or equal to the bulk lipid, not the
1179 effect of lengths exceeding it, and is only known for one bulk lipid
1180 component (DOPC). We assume then, that the increase is in relationship
1181 to the average vesicle, and that lipids with larger acyl chain length
1182 will also show an increase in their dissociation constant.
1185 l_\mathrm{b} = 3.2^{\left|\left<l_\mathrm{vesicle}\right>-l_\mathrm{monomer}\right|}
1186 \label{eq:length_backward}
1189 The most common $\left<l_\mathrm{vesicle}\right>$ is around $17.75$,
1190 which leads to a range of $\Delta \Delta G^\ddagger$ from
1191 $\Sexpr{format(digits=3,to.kcal(3.2^abs(12-17.75)))}
1192 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with length 12
1193 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$
1194 for monomers with length near 18 to
1195 $\Sexpr{format(digits=3,to.kcal(3.2^abs(24-17.75)))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
1196 for monomers with length 24. The full range of possible values of
1197 $l_\mathrm{b}$ are shown in \cref{fig:lb_graph}.
1199 % (for methods? From McLean84LIB: The activation free energies and free
1200 % energies of transfer from self-micelles to water increase by 2.2 and
1201 % 2.1 kJ mol-' per methylene group, respectively. \RZ{see if we can use it
1202 % to justify arranging our changed in activating energy around 1
1207 <<lb_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1208 trellis.device(new=F,color=TRUE)
1209 trellis.par.set(list(axis.line =list(col="transparent")))
1210 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1211 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1212 grid <- expand.grid(x=seq(12,24,length.out=20),
1213 y=seq(12,24,length.out=20))
1214 grid$z <- 3.2^(abs(grid$x-grid$y))/1e6
1215 print(wireframe(z~x*y,grid,cuts=50,
1218 scales=list(arrows=F,col=1,distance=0.75),
1219 xlab=list(label=expression(paste(symbol("\341"),italic(l[vesicle]),symbol("\361"))),rot=30),
1220 ylab=list(label=expression(italic(l[monomer])),rot=-35),
1221 zlab=list(label=expression(italic(l[b])%*%10^-6),rot=93)),newpage=FALSE)
1224 just=c("left","top"),
1225 gp=gpar(fontsize=48),
1226 y=unit(1,"npc")-unit(0.25,"lines"),
1227 x=unit(0,"npc")+unit(0.25,"lines"))
1229 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1230 grid <- expand.grid(x=seq(12,24,length.out=20),
1231 y=seq(12,24,length.out=20))
1232 grid$z <- to.kcal(3.2^(abs(grid$x-grid$y)))
1233 print(wireframe(z~x*y,grid,cuts=50,
1236 scales=list(arrows=F,col=1,distance=0.75),
1237 xlab=list(label=expression(paste(symbol("\341"),italic(l[vesicle]),symbol("\361"))),rot=30),
1238 ylab=list(label=expression(italic(l[monomer])),rot=-35),
1239 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1242 just=c("left","top"),
1243 gp=gpar(fontsize=48),
1244 y=unit(1,"npc")-unit(0.25,"lines"),
1245 x=unit(0,"npc")+unit(0.25,"lines"))
1248 \caption{Change in length backwards ($l_\mathrm{b}\times 10^{-6}$) (A) or activation energy
1249 in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the length of the
1250 monomer leaving ($l_\mathrm{monomer}$) and the average length of
1251 lipids in the vesicle ($\left<l_\mathrm{vesicle}\right>$).}
1252 \label{fig:lb_graph}
1256 \subsubsection{Complex Formation Backward}
1258 Complex formation ($CF1$) describes the interaction between \ac{CHOL} and
1259 \ac{PC} or \ac{SM}, where \ac{PC} or \ac{SM} protects the hydroxyl group of \ac{CHOL} from
1260 interactions with water %
1261 \citep{Huang1999:chol_solubility_model,
1262 Huang1999:cholesterol_solubility, McConnell2006, McConnell2003}%
1263 . \ac{PC} ($CF1=2$) can interact with two \ac{CHOL}, and \ac{SM} ($CF1=3$) with three
1264 \ac{CHOL} ($CF1=-1$). If the average of $CF1$ is positive (excess of \ac{SM} and
1265 \ac{PC} with regards to complex formation), components with negative $CF1$
1266 (\ac{CHOL}) will be retained. If average $CF1$ is negative, components with
1267 positive $CF1$ are retained. An equation which has this property is
1268 $CF1_\mathrm{b}=a^{\left<CF1_\mathrm{vesicle}\right>
1269 CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{vesicle}\right>
1270 CF1_\mathrm{monomer}\right|}$, where difference of the exponent is
1271 zero if the average $CF1$ and the $CF1$ of the component have the same
1272 sign, or double the product if the signs are different. Based on
1273 activation energy considerations, we took an ad hoc base for $a$ as
1278 CF1_\mathrm{b}=1.5^{\left<CF1_\mathrm{vesicle}\right> CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{vesicle}\right> CF1_\mathrm{monomer}\right|}
1279 \label{eq:complex_formation_backward}
1282 The most common $\left<CF1_\mathrm{vesicle}\right>$ is around $0.925$,
1283 which leads to a range of $\Delta \Delta G^\ddagger$ from
1284 $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*-1-abs(0.925*-1))))}
1285 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
1287 $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*2-abs(0.925*2))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
1288 for monomers with complex formation $2$ to
1289 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
1290 formation $0$. The full range of possible values for $CF1_\mathrm{b}$ are
1291 depicted in \cref{fig:cf1b_graph}.
1296 <<cf1b_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1297 trellis.device(new=F,color=TRUE)
1298 trellis.par.set(list(axis.line =list(col="transparent")))
1299 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1300 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1301 grid <- expand.grid(x=seq(-1,3,length.out=20),
1302 y=seq(-1,3,length.out=20))
1303 grid$z <- 1.5^(grid$x*grid$y-abs(grid$x*grid$y))
1304 print(wireframe(z~x*y,grid,cuts=50,
1307 scales=list(arrows=F,col=1,distance=0.75),
1308 xlab=list(label=expression(paste(symbol("\341"),italic(CF1[vesicle]),symbol("\361"))),rot=30),
1309 ylab=list(label=expression(italic(CF1[monomer])),rot=-35),
1310 zlab=list(label=expression(italic(CF1[b])),rot=93)),newpage=FALSE)
1313 just=c("left","top"),
1314 gp=gpar(fontsize=48),
1315 y=unit(1,"npc")-unit(0.25,"lines"),
1316 x=unit(0,"npc")+unit(0.25,"lines"))
1318 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1319 grid <- expand.grid(x=seq(-1,3,length.out=20),
1320 y=seq(-1,3,length.out=20))
1321 grid$z <- to.kcal(1.5^(grid$x*grid$y-abs(grid$x*grid$y)))
1322 print(wireframe(z~x*y,grid,cuts=50,
1325 scales=list(arrows=F,col=1,distance=0.75),
1326 xlab=list(label=expression(paste(symbol("\341"),italic(CF1[vesicle]),symbol("\361"))),rot=30),
1327 ylab=list(label=expression(italic(CF1[monomer])),rot=-35),
1328 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1331 just=c("left","top"),
1332 gp=gpar(fontsize=48),
1333 y=unit(1,"npc")-unit(0.25,"lines"),
1334 x=unit(0,"npc")+unit(0.25,"lines"))
1337 \caption{Change in complex formation 1 backwards ($CF1_\mathrm{b}$) (A) or
1338 activation energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus
1339 the complex formation of the monomer leaving
1340 ($CF1_\mathrm{monomer}$) and the average complex formation of the
1341 vesicle ($\left<CF1_\mathrm{vesicle}\right>$).}
1342 \label{fig:cf1b_graph}
1345 \section{Simulation Methodology}
1347 \subsection{Overall Architecture}
1349 The simulations are currently run using a program written in perl with
1350 various modules to handle the subsidiary parts. It produces output for
1351 each generation, including the step immediately preceding and
1352 immediately following a vesicle split (and optionally, each step) that
1353 is written to a state file which contains the state of the vesicle,
1354 environment, kinetic parameters, program invocation options, time, and
1355 various other parameters necessary to recreate the state vector at a
1356 given time. This output file is then read by a separate program in
1357 perl to produce different output which is generated out-of-band;
1358 output which includes graphs and statistical analysis is performed
1359 using R~\citep{R:maincite} (and various grid graphics modules) called
1360 from the perl program.
1362 The separation of simulation and output generation allows refining
1363 output, and simulations performed on different versions of the
1364 underlying code to be compared using the same analysis methods and
1365 code. It also allows later simulations to be restarted from a specific
1366 generation, utilizing the same environment.
1368 Random variables of different distributions are calculated using
1369 Math::Random (0.71), which is seeded for each run using entropy from
1370 the linux kernel's urandom device.
1372 Code is available upon request.
1374 \subsection{Environment Creation}
1376 \subsubsection{Components}
1377 The environment contains concentrations of components. In the current
1378 set of simulations, there are \Sexpr{1+4*5*7} different possible
1379 components, consisting of \ac{PC}, PE, \ac{PS}, \ac{SM}, and \ac{CHOL}, with all lipids
1380 except for \ac{CHOL} having 5 possible unsaturations ranging from 0 to 4,
1381 and 7 lengths from $12,14,...,22$ ($7\cdot
1382 5\cdot4+1=\Sexpr{1+4*5*7}$). In cases where the environment has less
1383 than the maximum number of components, the components are selected in
1384 order without replacement from a randomly shuffled deck of components
1385 (with the exception of \ac{CHOL}) represented once until the desired number
1386 of unique components is obtained. \ac{CHOL} is overrepresented
1387 $\Sexpr{5*7}$ times to be at the level of other lipid types, ensuring
1388 that the probability of \ac{CHOL} being absent in the environment is the
1389 same as the probability of any one of the other lipid types (\ac{PS}, PE,
1390 etc.) being absent. This reduces the number of simulations with a
1391 small number of components which are devoid of \ac{CHOL}.
1393 \subsubsection{Concentrations}
1394 Once the components of the environment have been selected, their
1395 concentrations are determined. In experiments where the environmental
1396 concentration is the same across all lipid components, the
1397 concentration is set to $10^{-10}$~M. In other cases, the
1398 environmental concentration is set to a random number from a gamma
1399 distribution with shape parameter 2 and an average of
1400 $10^{-10}$~M. The base concentration ($10^{-10}$~M)
1401 can also be altered at the initialization of the experiment to
1402 specific values for specific lipid types or components.
1404 The environment is a volume which is the maximum number of vesicles
1405 from a single simulation (4096) times the size of the vesicle
1406 (usually 10000) divided by Avagadro's number divided by the total
1407 environmental lipid concentration, or usually
1408 \Sexpr{4096*10000/6.022E23/141E-10}~L.
1410 \subsection{Initial Vesicle Creation}
1412 Initial vesicles are seeded by selecting lipid molecules from the
1413 environment until the vesicle reaches a specific starting size. The
1414 vesicle starting size has gamma distribution with shape parameter 2
1415 and a mean of the per-simulation specified starting size, with a
1416 minimum of 5 lipid molecules, or can be specified to have a precise
1417 number of molecules. Lipid molecules are then selected to be added to
1418 the lipid membrane according to four different methods. In the
1419 constant method, molecules are added in direct proportion to their
1420 concentration in the environment. The uniform method adds molecules in
1421 proportion to their concentration in the environment scaled by a
1422 uniform random value, whereas the random method adds molecules in
1423 proportion to a uniform random value. The final method is a binomial
1424 method, which adjust the probability of adding a molecule of a
1425 specific component by the concentration of that component, and then
1426 adds the molecules one by one to the membrane. This final method is
1427 also used in the first three methods to add any missing molecules to
1428 the starting vesicle which were unallocated due to the requirement to
1429 add integer numbers of molecules. (For example, if a vesicle was to
1430 have 10 molecules split evenly between three lipid types, the 10th
1431 molecule would be assigned by randomly choosing between the three
1432 lipid types, weighted by their concentration in the environment.)
1434 \subsection{Simulation Step}
1436 Once the environment has been created and the initial vesicle has been
1437 formed, molecules join and leave the vesicle based on the kinetic
1438 parameters and state equation discussed above until the vesicle splits
1439 forming two progeny vesicles. The program then follows both vesicles
1440 and their progeny until the simulation reaches either the maximum
1441 number of vesicles (usually 4096), or the maximum simulation time
1444 \subsubsection{Calculation of Vesicle Properties}
1446 \label{sec:ves_prop}
1447 $S_\mathrm{vesicle}$ is the surface area of the vesicle, and is the sum of
1448 the surface area of all of the individual lipid components; each lipid
1449 type has a different surface area; we assume that the lipid packing is
1450 optimal, and there is no empty space.
1452 For the other vesicle properties (length, unsaturation, charge, and
1453 curvature), we calculate the mean, the standard deviation, the mean of
1454 the log, the mean of the absolute value of the log, the standard
1455 deviation of the log, and the standard deviation of the absolute value
1456 of the log. All cases, when we take the log, we take the log of the
1457 absolute value, and map $\log(0)$ to $0$. For the purpose of plotting
1458 vesicle properties, we only plot the mean of the property.
1460 \subsubsection{Joining and Leaving of Lipid Molecules}
1462 Determining the number of molecules to add to the lipid membrane
1463 ($n_i$) requires knowing $k_{\mathrm{f}i_\mathrm{adj}}$, the surface area of the
1464 vesicle $S_\mathrm{vesicle}$ (see \cref{sec:ves_prop}), the time interval
1465 $dt$ during which lipids are added, the base $k_{\mathrm{f}i}$, and the
1466 concentration of the monomer in the environment
1467 $\left[C_{i\mathrm{monomer}}\right]$ (see \cref{eq:state}).
1468 $k_{\mathrm{f}i\mathrm{adj}}$ is calculated (see \cref{eq:kf_adj}) based on the
1469 vesicle properties and their hypothesized effect on the rate (in as
1470 many cases as possible, experimentally based)
1471 (see \cref{sec:kinetic_adjustments} for details). $dt$ can be varied
1472 (see \cref{sec:step_duration}), but for a given step is constant. This
1473 leads to the following:
1475 $n_i = k_{\mathrm{f}i}k_{\mathrm{f}i_\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]S_\mathrm{vesicle}\mathrm{N_A}dt$
1477 In the cases where $n_i > 1$, the integer number of molecules is
1478 added. Fractional $n_i$ or the fractional remainder after the addition
1479 of the integer molecules are the probability of adding a specific
1480 molecule, and are compared to a uniformly distributed random value
1481 between 0 and 1. If the random value is less than or equal to the
1482 fractional part of $n_i$, an additional molecule is added.
1484 Molecules leaving the vesicle are handled in a similar manner, with
1486 $n_i = k_{\mathrm{b}i}k_{\mathrm{b}i_\mathrm{adj}}C_{i_\mathrm{vesicle}}\mathrm{N_A}dt$.
1488 Before addition or removal of molecules, the properties of the vesicle
1489 are calculated; this is done so that molecules can be considered to be
1490 added and removed at the same instant, even though additions are
1491 handled first programmatically. This also avoids cases where a removal
1492 would have resulted in a negative number of molecules for a particular
1495 \subsubsection{Step duration}
1496 \label{sec:step_duration}
1498 $dt$ is the time taken for each step of joining, leaving, and checking
1499 split. In the current implementation, it starts with a value of
1500 $10^{-6}$~$\mathrm{s}$ but this is modified in between each step if the
1501 number of molecules joining or leaving is too large or small. If more
1502 than half of the starting vesicle size molecules join or leave in a
1503 single step, $dt$ is reduced by half. If less than the starting
1504 vesicle size molecules join or leave in 100 steps, $dt$ is doubled.
1505 This is necessary to curtail run times and to automatically adjust to
1506 different experimental runs.
1508 \subsubsection{Vesicle splitting}
1510 If a vesicle has grown to a size which is more than double the
1511 starting vesicle size, the vesicle splits. Molecules are assigned to
1512 the progeny vesicles at random, with each progeny vesicle having an
1513 equal chance of getting a single molecule. The number of molecules to
1514 assign to the first vesicle has binomial distribution with a
1515 probability of an event in each trial of 0.5 and a number of trials
1516 equal to the number of molecules.
1520 The environment, initial vesicle, and the state of the vesicle
1521 immediately before and immediately after splitting are stored
1522 to produce later output.
1524 % \section{Analyzing output}
1526 % The analysis of output is handled by a separate perl program which
1527 % shares many common modules with the simulation program. Current output
1528 % includes simulation progress, summary tables, summary statistics, and
1531 % \subsection{PCA plots}
1533 % Two major groups of axes that contribute to vesicle variation between
1534 % subsequent generations are the components and properties of vesicles.
1535 % Each component in a vesicle is an axis on its own; it can be measured
1536 % either as an absolute number of molecules in each component, or the
1537 % fraction of molecules of that component over the total number of
1538 % molecules; the second approach is often more convenient, as it allows
1539 % vesicles with different number of molecules to be directly compared
1540 % (though it hides any effect of vesicle size).
1542 % In order to visualize the variation between and transition of
1543 % subsequent generations of vesicles from their initial state in the
1544 % simulation, to their final state at the termination of the simulation,
1545 % we plot the projection of the generations onto the two principle PCA
1546 % axes (and alternatively, any pairing of the axes).
1548 % \subsection{Carpet plots}
1550 % Carpet plots show the distance/similarity between the composition or
1551 % properties of all generations in a single run. The difference between
1552 % each group of vesicle can be calculated using Euclidean distance
1553 % (properties and compositions) or H similarity (composition only). We
1554 % must use Euclidean distance for properties because the H distance
1555 % metric is invalid when comparing negative vector elements to positive
1558 % In addition to showing distance or similarity, carpet plots also
1559 % depict propersomes and composomes as square boxes on the diagonals and
1560 % propertypes and compotypes as rectangles off the diagonals, each
1561 % propertype or compotype with a distinct color.
1563 % \subsection{Propersomes, propertypes, composomes and compotypes}
1565 % A generation is considered to be a propersome if it is less than $0.6$
1566 % units (by Euclidean distance of normalized properties) away from the
1567 % generation immediately following and preceding. Likewise, a generation
1568 % is in a composome if its H similarity is more than $0.9$ (by the
1569 % normalized H metric) from the preceding and following generations.
1570 % Propersomes and composomes are then assigned to propertypes and
1571 % compotypes using Paritioning Around Medioids
1572 % (PAM). All values of $k$ between 2 and 15
1573 % (or the number of propersomes and composomes, if that is less than 15)
1574 % are tried, and the clustering with the smallest
1575 % silhouette~\citep{Rousseeuw1987:silhouettes} is chosen as the ideal
1576 % clustering~\citep{Shenhav2005:pgard}.
1578 \section*{Formalism}
1580 The most current revision of this formalism is available at
1581 \url{https://git.donarmstrong.com/ool/lipid_simulation_formalism.git}.
1582 This document is \gitMarkPref • \gitMark.
1584 %\bibliographystyle{unsrtnat}
1585 %\bibliography{references.bib}