1 \documentclass[english,12pt]{article}
3 %\usepackage[pdftex]{graphicx}
5 \usepackage[bf]{caption2}
13 \usepackage[light,all]{draftcopy}
19 \usepackage[noblocks]{authblk}
20 \usepackage[hyperfigures,bookmarks,colorlinks,citecolor=black,filecolor=black,linkcolor=black,urlcolor=black]{hyperref}
21 \usepackage[sectionbib,sort&compress,numbers]{natbib}
22 \usepackage[nomargin,inline,draft]{fixme}
23 \usepackage[x11names,svgnames]{xcolor}
25 \newenvironment{narrow}[2]{%
27 \setlength{\topsep}{0pt}%
28 \setlength{\leftmargin}{#1}%
29 \setlength{\rightmargin}{#2}%
30 \setlength{\listparindent}{\parindent}%
31 \setlength{\itemindent}{\parindent}%
32 \setlength{\parsep}{\parskip}}%
34 \newenvironment{paperquote}{%
39 \renewcommand{\textfraction}{0.15}
40 \renewcommand{\topfraction}{0.85}
41 \renewcommand{\bottomfraction}{0.65}
42 \renewcommand{\floatpagefraction}{0.60}
43 %\renewcommand{\baselinestretch}{1.8}
44 \newenvironment{enumerate*}%
46 \setlength{\itemsep}{0pt}%
48 \setlength{\parskip}{0pt}}%
50 \newenvironment{itemize*}%
52 \setlength{\itemsep}{0pt}%
53 \setlength{\parskip}{0pt}}%
56 \newcommand{\DLA}[1]{\textcolor{red}{\fxnote{DLA: #1}}}
58 \newcommand{\RZ}[1]{\textcolor{red}{\fxnote{RZ: #1}}}
60 \newcommand{\OM}[1]{\textcolor{red}{\fxnote{OM: #1}}}
70 \title{Supplemental Material}
72 % rubber: module bibtex
73 % rubber: module natbib
78 <<load.libraries,echo=FALSE,results="hide",warning=FALSE,message=FALSE,error=FALSE,cache=FALSE>>=
79 opts_chunk$set(dev="CairoPDF",out.width="\\columnwidth",out.height="0.7\\textheight",out.extra="keepaspectratio")
80 opts_chunk$set(cache=TRUE, autodep=TRUE)
81 options(scipen = -2, digits = 1)
86 to.latex <- function(x){
87 gsub("\\\\","\\\\\\\\",latexSN(x))
90 to.kcal <- function(k,temp=300) {
92 return(-gasconst*temp*log(k)/1000)
96 % \setcounter{figure}{0} \setcounter{table}{0}
98 \renewcommand{\thefigure}{S\@arabic\c@figure}
99 \renewcommand{\thetable}{S\@arabic\c@table}
102 \section{Competition Implementation}
103 \subsection{Implementation changes}
106 \item settable maximum number of vesicles to track (default $10^4$)
107 \item start with 1~L ($10^{-3}$~m$^3$) cube
108 \item if at any point the number of vesicles exceeds the maximum
109 number, chop the volume and environment molecule number into tenths,
110 randomly select one tenth of the vesicles, and continue tracking.
111 \item generations will be counted per vesicle, and each progeny
112 vesicle will have a generation number one greater than the parental
114 \item 100 generations can result in as many as $2^{100}$
115 ($\Sexpr{to.latex(format(digits=3,2^100))}$) vesicles or as few as
117 \item Environment will use a specific number of each component instead
118 of a constant concentration; as the number may be larger than
119 \texttt{long long} ($2^{64}$), we use libgmp to handle an arbitrary
120 precision number of components
123 \subsection{Infrastructure changes}
125 \item Rewrite core bits in C
126 \item Use libgmp for handling large ints
127 \item Use openmpi to split the calculations out over multiple
128 machines/processors and allow deploying to larger
129 clusters/supercomputers
134 \section{State Equation}
135 % double check this with the bits in the paper
137 For a system with monomers $(_\mathrm{monomer})$ and a vesicle
138 $(_\mathrm{vesicle})$, the change in composition of the $i$th component of
139 a lipid vesicle per change in time ($d \left[C_{i_\mathrm{vesicle}}\right]/dt$)
140 can be described by a modification of the basic mass action law:
143 \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] -
144 k_{\mathrm{b}i}k_{\mathrm{b}i\mathrm{adj}}\left[C_{i_\mathrm{vesicle}}\right]
148 The base forward kinetic parameter for the $i$th component is $k_{\mathrm{f}i}$
149 and is dependent on the particular lipid type (PC, PS, SM, etc.). The
150 forward adjustment parameter, $k_{\mathrm{f}i\mathrm{adj}}$, is based on the
151 properties of the vesicle and the specific component (type, length,
152 unsaturation, etc.) (see Equation~\ref{eq:kf_adj}, and
153 Section~\ref{sec:kinetic_adjustments}).
154 $\left[C_{i_\mathrm{monomer}}\right]$ is the molar concentration of
155 monomer of the $i$th component. $\left[S_\mathrm{vesicle}\right]$ is the surface
156 area of the vesicle per volume. The base backwards kinetic parameter
157 for the $i$th component is $k_{\mathrm{b}i}$ and its adjustment parameter
158 $k_{\mathrm{b}i\mathrm{adj}}$ (see Equation~\ref{eq:kb_adj}, and
159 Section~\ref{sec:kinetic_adjustments}).
160 $\left[C_{i_\mathrm{vesicle}}\right]$ is the molar concentration of
161 the $i$th component in the vesicle.
163 \subsection{Per-Lipid Kinetic Parameters}
165 <<echo=FALSE,results='hide'>>=
166 kf.prime <- c(3.7e6,3.7e6,5.1e7,3.7e6,2.3e6)
167 kf <- (as.numeric(kf.prime)*10^-3)/(63e-20*6.022e23)
170 Each of the 5 lipid types has different kinetic parameters; where
171 available, these were taken from literature.
175 % \begin{tabular}{c c c c c c c c}
177 % Type & $k_\mathrm{f}$ $\left(\frac{\mathrm{m}}{\mathrm{s}}\right)$ & $k'_\mathrm{f}$ $\left(\frac{1}{\mathrm{M} \mathrm{s}}\right)$ & $k_\mathrm{b}$ $\left(\mathrm{s}^{-1}\right)$ & Area $\left({\AA}^2\right)$ & Charge & $\mathrm{CF}1$ & Curvature \\
179 % PC & $\Sexpr{to.latex(format(digits=3,scientific=TRUE,kf[1]))}$ & $3.7 \times 10^6$ & $2 \times 10^{-5}$ & 63 & 0 & 2 & 0.8 \\
180 % PS & $\Sexpr{to.latex(format(digits=3,scientific=TRUE,kf[2]))}$ & $3.7 \times 10^6$ & $1.25\times 10^{-5}$ & 54 & -1 & 0 & 1 \\
181 % CHOL & $\Sexpr{to.latex(format(digits=3,scientific=TRUE,kf[3]))}$ & $5.1 \times 10^7$ & $2.8 \times 10^{-4}$ & 38 & 0 & -1 & 1.21 \\
182 % SM & $\Sexpr{to.latex(format(digits=3,scientific=TRUE,kf[4]))}$ & $3.7 \times 10^6$ & $3.1 \times 10^{-3}$ & 61 & 0 & 3 & 0.8 \\
183 % PE & $\Sexpr{to.latex(format(digits=3,scientific=TRUE,kf[5]))}$ & $2.3 \times 10^6$ & $1 \times 10^{-5}$ & 55 & 0 & 0 & 1.33 \\
186 % \caption{Kinetic parameters and molecular properties of lipid types}
187 % \label{tab:kinetic_parameters_lipid_types}
190 %%% \DLA{I think we may just reduce these three sections; area, $k_\mathrm{f}$
191 %%% and $k_\mathrm{b}$ to Table~\ref{tab:kinetic_parameters_lipid_types} with
194 %%% \RZ{Yes, but we also have to have then as comments the numbers that
195 %%% are not supported by refs}
197 \subsubsection{$k_\mathrm{f}$ for lipid types}
198 $k_{\mathrm{f}_\mathrm{PC}}$ was measured by
199 \citet{Nichols1985:phospholipid_monomer_vesicle_thermodynamics} and
200 was found to be $3.7\times 10^6$~$\frac{1}{\mathrm{M}\mathrm{s}}$ by the
201 partitioning of P-C$_6$-NBD-PC between DOPC vesicles and water. As
202 similar references are not available for SM or PS, we assume that they have
203 the same $k_\mathrm{f}$. For CHOL, no direct measurement of $k_\mathrm{f}$ is available,
204 however, \citet{Estronca2007:dhe_kinetics} measured the transfer of
205 DHE from BSA to LUV, and found a $k_\mathrm{f}$ of $5.1\times 10^7$~%
206 $\frac{1}{\mathrm{M} \mathrm{s}}$. We assume that this value is close
207 to that of CHOL, and use it for $k_{\mathrm{f}_\mathrm{CHOL}}$. In the case of
208 PE, \citet{Abreu2004:kinetics_ld_lo} measured the association of
209 NBD-DMPE with POPC LUV found a value for $k_\mathrm{f}$ of $2.3 \times 10^{6}$~%
210 $\frac{1}{\mathrm{M} \mathrm{s}}$. These three authors used a slightly
211 different kinetic formalism ($\frac{d\left[A_v\right]}{dt} =
212 k'_\mathrm{f}[A_m][L_v] - k_\mathrm{b}[A_v]$), so we correct their values of $k_\mathrm{f}$ by
213 multiplying by the surface area of a mole of lipids.
215 \subsubsection{$k_\mathrm{b}$ for lipid types}
217 \citet{Wimley1990:dmpc_exchange} measured the half time for the
218 exchange of [$^3$H]DMPC between LUV at 30\textdegree C, and found it
219 to be 9.6 hr. As this is a first order reaction, and the primary
220 limiting step in exchange is lipid desorption, $k_\mathrm{b}$ for DMPC is
221 $k_{\mathrm{b}_\mathrm{PC}}=\frac{\log 2}{9.6 \times 60 \times 60} \approx
222 \Sexpr{log(2)/(9.6*60*60)}
223 \mathrm{s}^{-1}$. We assume that $k_\mathrm{b}$ for SM is the same as for PC.
224 To estimate the $k_\mathrm{b}$ of PE and PS, we used the data from
225 \citet{Nichols1982:ret_amphiphile_transfer} who measured the rate of
226 exchange of the fluorescent label C$_6$-NBD attached to different
227 lipid species. Although the values of kb are different for the labeled
228 and unlabeled lipids, we assume that the ratios of the kinetics
229 constants for the lipid types are the same. Furthermore we assume that
230 PG behaves similarly to PS. Thus, we can determine the $k_\mathrm{b}$ of PE and
231 PS from the already known $k_\mathrm{b}$ of PC. For C$_6$-NBD labeled PC,
232 \citet{Nichols1982:ret_amphiphile_transfer} obtained a $k_\mathrm{b}$ of
233 $0.89$~$\mathrm{min}^{-1}$, PE of $0.45$~$\mathrm{min}^{-1}$ and PG of
234 $0.55$~$\mathrm{min}^{-1}$. Assuming the ratios of the rate of
235 exchange are the same for unlabeled lipids and labeled lipids, we can
236 determine the $k_\mathrm{b}$ of PE and PS from the already known $k_\mathrm{b}$ of
237 PC~\citep{Wimley1990:dmpc_exchange}. Calculating the ratio, this leads
238 us to $k_{\mathrm{b}_\mathrm{PE}} =
239 \frac{k_{\mathrm{b}_\mathrm{PC}}\times\mathrm{PE}}{\mathrm{PC}} \approx
240 \frac{2\times 10^{-5}\,\mathrm{s}^{-1} \times
241 0.45\,\mathrm{min}^{-1}}{0.89\,\mathrm{min}^{-1}} \approx
242 \Sexpr{log(2)/(9.6*60*60)*0.45/0.89}$~$\mathrm{s}^{-1}$
243 and likewise, $k_{\mathrm{b}_\mathrm{PS}}\approx
244 \Sexpr{log(2)/(9.6*60*60)*0.55/0.89}$~$\mathrm{s}^{-1}$.
245 The $k_\mathrm{b}$ of SM was determined using the work of
246 \citet{Bai1997:lipid_movementbodipy}, who measured spontaneous
247 transfer of C$_5$-DMB-SM and C$_5$-DMB-PC from donor and acceptor
248 vesicles, finding $3.4\times10^{-2}$~$\mathrm{s}^{-1}$ and
249 $2.2\times10^{-3}$~$\mathrm{s}^{-1}$ respectively; using the ratio of
250 $k_\mathrm{b}$ of C$_5$-DMB-SM to the $k_\mathrm{b}$ of C$_5$-DMB-PC times the $k_\mathrm{b}$ of
251 PC ($\frac{3.4 \times 10^{-2}\mathrm{s}^{-1}}{2.2 \times
252 10^{-3}\mathrm{s}^{-1}}\approx
253 \Sexpr{log(2)/(9.6*60*60)}\mathrm{s}^{-1}$),
254 we obtain $k_{\mathrm{b}_\mathrm{SM}} \approx
255 \Sexpr{log(2)/(9.6*60*60)*3.4e-2/2.2e-3}$.
256 In the case of CHOL, \citet{Jones1990:lipid_transfer} measured the
257 $t_{1/2}$ of [$^3$H] transfer from POPC vesicles and found it to be 41
258 minutes, leading to a $k_{\mathrm{b}_\mathrm{CHOL}} = \frac{\log 2}{41\times
260 \Sexpr{log(2)/(41*60)}$~$\mathrm{s}^{-1}$.
262 \subsubsection{Headgroup Surface Area for lipid types}
264 % Smondyrev, A. M., and M. L. Berkowitz. 1999b. United atom force
265 % field for phospholipid membranes: constant pressure molecular
266 % dynamics simulation of dipalmitoylphosphatidicholine/water system.
267 % J. Comput. Chem. 50:531–545.
269 Different lipids have different headgroup surface areas, which contributes to
270 $\left[S_\mathrm{vesicle}\right]$. \citet{Smaby1997:pc_area_with_chol}
271 measured the surface area of POPC with a Langmuir film balance, and
272 found it to be 63~Ã…$^2$ at $30$~$\frac{\mathrm{mN}}{\mathrm{m}}$.
273 Molecular dynamic simulations found an area of 54 Ã…$^2$ for
274 DPPS\citep{Cascales1996:mds_dpps_area,Pandit2002:mds_dpps}, which is
275 in agreement with the experimental value of 56~Ã…$^2$ found using a
276 Langmuir balance by \citet{Demel1987:ps_area}.
277 \citet{Shaikh2002:pe_phase_sm_area} measured the area of SM using a
278 Langmuir film balance, and found it to be 61~Ã…$^2$. Using $^2$H NMR,
279 \citet{Thurmond1991:area_of_pc_pe_2hnmr} found the area of
280 DPPE-d$_{62}$ to be 55.4 Ã…$^2$. \citet{Robinson1995:mds_chol_area}
281 found an area for CHOL of 38~Ã…$^2$ using molecular dynamic
284 % robinson's chol area is kind of crappy; they did it using MDS, but
285 % vaguely; they don't indicate what it was at the end, and they only
286 % had 2 molecules of CHOL.
290 % PC: 63 A (Smaby97rd) (Pandit (?)
291 % PS: 54 A (Pandit02LIB) (Cascales 1996; J. Chem. Phys 104:2713, Mukhopadhyay2004)
292 % CHOL: 38 A (Robinson 1995) (Lundberg 1982)
293 % SM: 51 A (Shaikh2002, but 61(?))
294 % PE: 55 A (Thurmond, 1991) (Pandit2002)
296 % Compare to results by Petrache2004 with DOPC of 72.5, DOPS of 65.3.
298 % The cross‑sectional area of cholesterol is ~37 ֵ2 (Lundberg, B. 1982. A
299 % surface film study of the lateral packing of phosphatidylcholine and
300 % cholesterol. Chem. Phys. Lipids. 31:23-32).
302 % From Pandit02LIB: the area per headgroup for PE is ~10‑20\% (Thurmond
303 % et al., 1991) smaller than the area per PC molecule. For DLPE, it is
304 % 50.6 Öµ2 (McIntosh and Simon, 1986; Damodaran et al., 1992). The
305 % estimated area per molecule for DPPE is ~55.4 Öµ2 (Thurmond et al.,
308 % From Pandit02LIB: (from their simulations) average area per headgroup
309 % for DPPS molecules is 53.75 ± 0.10 ֵ2. The values inferred from the
310 % experiments are in the region between 45 and 55 Öµ2 (Cevc et al., 1981.
311 % Titration of the phase transition of phosphatidylserine bilayer
312 % membranes. Effect of pH, surface electrostatics, ion binding, and
313 % head‑group hydration. Biochemistry. 20:4955‑4965; | Demel et al.,
314 % 1987. Monolayer characteristics and thermal behavior of natural and
315 % synthetic phosphatidylserines. Biochemistry. 26:8659‑8665). Note that
316 % the average area per headgroup for DPPC bilayer obtained from
317 % simulations and measurements is ~62 Öµ A2. One would expect that the
318 % area per headgroup in case of DPPS molecules will be larger than the
319 % area per DPPC molecule, because the headgroups of DPPS are charged and
320 % therefore repel each other. Contrary to this expectation, the area per
321 % headgroup in DPPS is ~13\% smaller than that of DPPC. López Cascales
322 % et al. (1996. Molecular dynamics simulation of a charged biological
323 % membrane. J. Chem. Phys. 104:2713‑2720.) proposed that this reduction
324 % in the area per headgroup is due to a strong intermolecular
325 % coordination between DPPS molecules. MD gave area per POPS of 55 A2 at
326 % 300K (Mukhopadhyay04LIB).
328 % However, by 2H NMR and X-ray, Petrache04LIB determined the area of
329 % DOPS at 30 C to be 65.3 A2, considerably smaller than DOPC, determined
330 % by the same group in another paper to be 72.5 Öµ A2
331 % (Tristram‑Nagle00LIB).
333 % From Shaikh02: At 35 mN/m of surface pressure, the areas/molecule for
334 % POPE, SM, and CHOL were found to be 65.2, 61.1, and 38.8 Ã…2.
336 % In Lonchin99: POPC area 72 A2 ((a) Cornell, B. A.; Middlehurst, J.;
337 % Separovic, F. Biochim. Biophys. Acta1980, 598, 405−410), oleic acid
338 % area 32 A2 ((a) Small, D. M. In Handbook of Lipid Research; Hanahan,
339 % D. J., Ed.; Plenum Press:Â New York, 1986; Vol. 4, Chapter 9.
341 % (b) Jain, M. K. Introduction to Biological Membranes, 2nd ed.; John
342 % Wiley \& Sons:Â New York, 1988; Chapter 3). Cistola, D. P.; Atkinson,
343 % D.; Hamilton, J. A.; Small, D. M. Biochemistry1986, 25, 2804−2812.
345 % From Marsh96: The values obtained for the headgroup area per molecule,
346 % fitted as free parameters in Fig. 6, are APE = 58 A2 and APC = 78 A2,
347 % for DOPE and DOPC, respectively. These values are within the range
348 % obtained for the corresponding areas per molecule of
349 % phosphatidylethanolamines and phosphatidylcholines in the fluid
350 % lamellar phase (see, e.g., Marsh, 1990), and correspond to values for
351 % the lipid packing parameter of V/Al = 1.08 and 1.34 for DOPC and DOPE,
352 % respectively. The former value is consistent with DOPC alone forming
353 % hydrated lamellar structures rather than nonlamellar structures. The
354 % latter value lies at the lower end of the packing parameters obtained
355 % for phosphatidylethanolamines.
357 % From Kumar91: The head group area chosen was 71.7 A2 for PC and
358 % lysoPC, 42 A2 for PE, and 19 A2 for Chol (5, 14). The value selected
359 % for PC is in agreement with the 71 A2 value determined for PCs with
360 % saturated acyl chains (14). The hydrocarbon volume and the length
361 % chosen for Chol are 400 A3 and 17.25 A, respectively (5).
365 % From Sampaio05: Besides this work and our own earlier report on the
366 % association of NBD-DMPE with lipid bilayers (Abreu et al., 2004), we
367 % are aware of only one other report in the literature (Nichols, 1985)
368 % in which all the kinetic constants of lipid-derived amphiphile
369 % association with lipid bilayer membranes were experimentally measured.
370 % {indeed; everything is k- !!!; rz}
372 % From McLean84LIB: Although it is difficult to measure cmc values for
373 % the sparingly soluble lipids used in this study, estimates for
374 % lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds,
375 % 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution;
376 % Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6
377 % X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X
378 % 10-8 M for DMPC was estimated from the linear relationship between ln
379 % cmc and the number of carbons in the PC acyl chain by using data for n
380 % = 7, 8, 10, and 16 [summarized in Tanford (1980)].
382 % From Toyota08: Recently, several research groups have reported
383 % self-reproducing systems of giant vesicles that undergo a series of
384 % sequential transformations: autonomous growth, self-division, and
385 % chemical reactions to produce membrane constituents within the giant
388 % Vesicle sizes of 25 nm for SUV and 150 nm for LUV were mentioned by
391 % From Lund-Katz88: Charged and neutral small unilamellar vesicles
392 % composed of either saturated PC, unsaturated PC, or SM had similar
393 % size distributions with diameters of 23 \& 2 nm.
395 % From Sampaio05LIB: The exchange of lipids and lipid derivatives
396 % between lipid bilayer vesicles has been studied for at least the last
397 % 30 years. Most of this work has examined the exchange of amphiphilic
398 % molecules between a donor and an acceptor population. The measured
399 % efflux rates were shown in almost all cases, not surprisingly, to be
400 % first order processes. In all of this work, the insertion rate has
401 % been assumed to be much faster than the efflux rate. Having measured
402 % both the insertion and desorption rate constants for amphiphile
403 % association with membranes, our results show that this assumption is
404 % valid. In several cases reported in the literature, the insertion rate
405 % constant was assumed, although never demonstrated, to be a
406 % diffusion-controlled process.
408 % (for methods? From McLean84LIB: The activation free energies and free
409 % energies of transfer from self-micelles to water increase by 2.2 and
410 % 2.1 kJ mol-' per methylene group, respectively. {see if we can use it
411 % to justify arranging our changed in activating energy around 1
414 % Jones90 give diameter of LUV as 100 nm, and of SUV as 20 nm; that
415 % would give the number of molecules per outer leaflet of a vesicle as
418 % Form Simard08: Parallel studies with SUV and LUV revealed that
419 % although membrane curvature does have a small effect on the absolute
420 % rates of FA transfer between vesicles, the ΔG of membrane desorption
421 % is unchanged, suggesting that the physical chemical properties which
422 % govern FA desorption are dependent on the dissociating molecule rather
423 % than on membrane curvature. However, disagreements on this fundamental
424 % issue continue (57, 61, 63, 64)
426 % (methods regarding the curvature effect: Kleinfeld93 showed that the
427 % transfer parameters of long-chain FFA between the lipid vesicles
428 % depend on vesicle curvature and composition. Transfer of stearic acid
429 % is much slower from LUV as compared to SUV).
431 % From McLean84: In a well-defined experimental system consisting of
432 % unilamellar lipid vesicles, in the absence of protein, the
433 % rate-limiting step for the overall exchange process is desorption
434 % (McLean \& Phillips, 1981). {thus I can take exchange data for the
435 % estimation of k- rz; 8/11/08}.
437 \subsubsection{Complex Formation 1}
439 \citet{Thomas1988:chol_transfer} found that SM significantly decreases
440 the rate of cholesterol transfer from PC, while PE and PS have no
441 effect at physiologically significant levels. We attribute this to the
442 formation of complexes between SM and CHOL, which retards the rate of
443 efflux of CHOL from the membrane. Similar complexes exist between PC
444 and CHOL, where it was shown that two CHOL molecules complex with one
445 PC~\citep{Huang1999:chol_solubility_model,
446 Huang1999:cholesterol_solubility,McConnell2006,McConnell2003}. It
447 has also been shown that SM binds more CHOL molecules than
448 PC~\citep{Katz1988:pl_packing_chol}. We assume that three CHOL complex
449 with one SM, leading to $\mathrm{CF}1$ values of $3$, $2$, and $-1$
450 for SM, PC, and CHOL, respectively.
452 \subsubsection{Curvature}
454 We used the formulation for curvature of
455 \citet{Israelachvili1975:amphiphile_self_assembly}, or $S=\frac{v}{l_c
456 a}$, where $S$ is the curvature, which ranges from $0$ to $\infty$,
457 $l_c$ is the critical length of the acyl chain, $v$ is the volume of
458 the lipid, and $a$ is the area of the head group.
459 \citet{Kumar1991:lipid_packing} found a value for $S$ of 1.33 for PE,
460 1.21 for CHOL, and 0.8 for diC$_{16}$ PC. We assume that PS has neutral
461 curvature (value of 1), and that SM has the same curvature as PC
462 (0.8). It should also be noted that in reality the curvature parameter
463 changes with length, but at longer chain lengths, is effectively
464 constant; in the current model, curvature is held constant for each
468 % 1. Israelachvili, J. N., Mitchell, D. J. & Ninham, B. W. (1976) J.
469 % Chem. Soc. Faraday Trans. 2 72, 1525-1568.
470 % 2. Israelachvili, J. N., Marcelja, S. & Horn, R. G. (1980) Q. Rev.
471 % Biophys. 13, 121-200.
472 % 3. Israelachvili, J. N. (1985) Intermolecular and Surface Forces
473 % (Academic, New York), pp. 229-271.
475 % \DLA{This section needs to be added still.}
477 % PE = 1.33 (Kumar91)
478 % CHOL = 1.21 (Kumar91)
480 % SM=0.8 (assumed by rz same as PC)
481 % PS=1 (no refs so far; should be close to unity; rz)
483 % From Epand94LIB: a considerable amount of heat is evolved per mole of
484 % bilayer stabilizer when such molecules are inserted into membranes
485 % which are under a high curvature strain. If this energy were coupled
486 % to a membrane event, such as the conformational change in a protein,
487 % it could be the driving force responsible for such processes.
489 % From Epand94LIB: considerable energy may be released upon the
490 % incorporation of certain molecules into membranes which have a low
491 % radius of spontaneous curvature. {discuss with DL in terms of
492 % catalytic activity; rz; 8/24/2010}.
494 % Kleinfeld93 showed that the transfer parameters of long-chain FFA
495 % between the lipid vesicles depend on vesicle curvature and
496 % composition. Transfer of stearic acid is much slower from LUV as
499 % Substitution of DPPC SUV for POPC SUV as the donor matrix revealed
500 % some differences in the energetics of transfer that are associated
501 % with the physical state of the donor particle. Whereas the free
502 % energies of transfer of oleic acid from DPPC SUV and POPC SUV are
503 % similar, the free energy of transfer of stearic acid from DPPC SUV is
504 % much higher than that observed from POPC.
506 % From Kumar91: S is given by S = V/al, where V is the hydrocarbon
507 % volume, a is the area of head group, and l is the critical length of
508 % the hydrocarbon chain (1-3). a, V, and I are all estimable or
509 % measurable (4), and the value of S can be calculated. The value of S
510 % determines the aggregate formed by lipids or any amphiphiles upon
511 % hydration. It has been shown that lipids aggregate to form spherical
512 % micelles (S < 1/3), nonspherical (cylindrical) micelles (1/3 < S <
513 % 1/2), bilayers (1/2 < S < 1), and reverse micelles or hexagonal (HII,)
514 % phases (S > 1). However, theoreticians caution that the above
515 % predicted limits set on the values of S are relatively insensitive to
516 % the exact values of V and a but are strongly dependent upon the choice
519 % From Kumar91: Cylindrical and wedge-shaped molecules have been shown
520 % to be essential for spontaneous vesiculation (6) as well as bilayer
521 % stabilization (7, 8).
523 % From Kumar91: It has also been shown experimentally that molecules
524 % having complementary molecular shapes can form bilayer structures.
525 % Inverted cone (wedge)- shaped lysoPC and cone-shaped Chol or
526 % unsaturated PE can interact stoichiometrically to form bilayer
527 % structures (9-12). {quote to justify our formula for complementary
530 \section{Kinetic adjustments}
531 \label{sec:kinetic_adjustments}
533 In the mass action equation used in the formalism, both the forward
534 and backwards kinetic constants ($k_\mathrm{f}$ and $k_\mathrm{b}$) are adjusted (by
535 $k_{\mathrm{f}i\mathrm{adj}}$ and $k_{\mathrm{b}i\mathrm{adj}}$) to reflect the
536 properties of the vesicle and the properties of the monomer entering
537 or exiting the vesicle.
539 \subsection{Forward Kinetic Adjustments}
541 The forward rate constant adjustment, $k_{\mathrm{f}i\mathrm{adj}}$ takes into
542 account unsaturation ($un_\mathrm{f}$), charge ($ch_\mathrm{f}$), curvature ($cu_\mathrm{f}$),
543 length ($l_\mathrm{f}$), and complex formation ($CF1_\mathrm{f}$), each of which is
544 modified depending on the specific component and the vesicle into
545 which the component is entering; the final $k_{\mathrm{f}i\mathrm{adj}}$ is the
546 product of the adjustments made separately for each property.
549 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}
553 \subsubsection{Unsaturation Forward}
555 In order for a lipid to be inserted into a membrane, a void has to be
556 formed for it to fill. Voids can be generated by the combination of
557 unsaturated and saturated lipids forming heterogeneous domains. Void
558 formation is increased when the unsaturation of lipids in the vesicle
559 is widely distributed; in other words, the insertion of lipids into
560 the membrane is greater when the standard deviation of the
561 unsaturation is larger %
562 %%% \RZ{May need to site (at least for us) works showing
563 %%% mismatch-dependent ``defects''}. %
564 Assuming that an increase in width of the distribution linearly
565 decreases the free energy of activation, the $un_\mathrm{f}$ parameter must
566 follow $a^{\mathrm{stdev}\left(un_\mathrm{vesicle}\right)}$ where $a > 1$,
567 so a convenient starting base for $a$ which leads to a reasonable
568 change in Eyring activation energy is $2$:
571 un_\mathrm{f} = 2^{\mathrm{stdev}\left(un_\mathrm{vesicle}\right)}
572 \label{eq:unsaturation_forward}
575 The mean $\mathrm{stdev}\left(un_\mathrm{vesicle}\right)$ in our
576 simulations is around $1.5$, which leads to a $\Delta \Delta
577 G^\ddagger$ of $\Sexpr{to.latex(format(digits=3,to.kcal(2^1.5)))}
578 \frac{\mathrm{kcal}}{\mathrm{mol}}$, and a total range of possible
579 values depicted in Figure~\ref{fig:unf_graph}.
581 % \RZ{Explain here, or even earlier that the formulas were ad hoc
582 % adjusted to correspond to ``reasonable'' changes in the Eyring
583 % activation energies.}
585 % It is not clear that the unsaturation of the inserted monomer will
586 % affect the rate of the insertion positively or negatively, so we do
587 % not include a term for it in this formalism.
591 <<unf_graph,echo=FALSE,results='hide',fig.width=10,fig.height=5,cache=FALSE>>=
592 layout(matrix(1:2,nrow=1,ncol=2))
593 curve(2^x,from=0,to=sd(c(0,4)),
594 xlab=expression(paste(stdev,group("(",italic(un[vesicle]),")"))),
595 ylab=expression(italic(un[f])))
596 vps <- baseViewports()
597 pushViewport(vps$figure)
599 just=c("left","top"),
600 gp=gpar(fontsize=48),
601 y=unit(1,"npc")-unit(0.25,"lines"),
602 x=unit(0,"npc")+unit(0.25,"lines"))
604 curve(to.kcal(2^x),from=0,to=sd(c(0,4)),
605 xlab=expression(paste(stdev,group("(",italic(un[vesicle]),")"))),
606 ylab="Activation Energy (kcal/mol)")
607 vps <- baseViewports()
608 pushViewport(vps$figure)
610 just=c("left","top"),
611 gp=gpar(fontsize=48),
612 y=unit(1,"npc")-unit(0.25,"lines"),
613 x=unit(0,"npc")+unit(0.25,"lines"))
616 \caption{Change in unsaturation forward ($un_\mathrm{f}$) (A) or activation
617 energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the
618 standard deviation of the unsaturation of the vesicle
619 ($\mathrm{stdev}\left(un_\mathrm{vesicle}\right)$).}
620 \label{fig:unf_graph}
623 \subsubsection{Charge Forward}
625 A charged lipid such as PS approaching a vesicle with an average
626 charge of the same sign will experience repulsion, whereas those with
627 different signs will experience attraction, the degree of which is
628 dependent upon the charge of the monomer and the average charge of the
629 vesicle. If either the vesicle or the monomer has no charge, there
630 should be no effect of charge upon the rate. This leads us to the
631 following equation, $a^{-\left<ch_v\right> ch_m}$ (which is similar to
632 the numerator of Coulomb's law), where $\left<ch_v\right>$ is the
633 average charge of the vesicle, and $ch_m$ is the charge of the
634 monomer. If either $\left<ch_v\right>$ or $ch_m$ is 0, the adjustment
635 parameter is 1 (no change), whereas it decreases if both are positive
636 or negative, as the product of two real numbers with the same sign is
637 always positive. The base for $a$ ($60$) was chosen ad hoc to
638 correspond to a maximum of about
639 $0.5\frac{\mathrm{kcal}}{\mathrm{mol}}$ change in activation energy
640 for the common case, resulting in the following equation:
644 ch_\mathrm{f} = 60^{-\left<{ch}_v\right> {ch}_m}
645 \label{eq:charge_forward}
648 The most common $\left<{ch}_v\right>$ is around $-0.165$, which leads
649 to a range of $\Delta \Delta G^\ddagger$ from
650 $\Sexpr{format(digits=3,to.kcal(60^(-.165*-1)))}
651 \frac{\mathrm{kcal}}{\mathrm{mol}}$ to
652 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$, and the total range of possible
653 values seen in Figure~\ref{fig:chf_graph}.
657 <<chf_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
658 trellis.device(new=F,color=TRUE)
659 trellis.par.set(list(axis.line =list(col="transparent")))
660 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
661 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
662 x <- seq(-1,0,length.out=20)
663 y <- seq(-1,0,length.out=20)
664 grid <- expand.grid(x=x,y=y)
665 grid$z <- as.vector(60^(-outer(x,y)))
666 print(wireframe(z~x*y,grid,cuts=50,
669 scales=list(arrows=F,col=1,distance=0.75),
670 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
671 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
672 zlab=list(label=expression(italic(ch[f])),rot=93)),newpage=FALSE)
675 just=c("left","top"),
676 gp=gpar(fontsize=48),
677 y=unit(1,"npc")-unit(0.25,"lines"),
678 x=unit(0,"npc")+unit(0.25,"lines"))
680 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
681 x <- seq(-1,0,length.out=20)
682 y <- seq(-1,0,length.out=20)
683 grid <- expand.grid(x=x,y=y)
684 grid$z <- as.vector(to.kcal(60^(-outer(x,y))))
685 print(wireframe(z~x*y,grid,cuts=50,
688 scales=list(arrows=F,col=1,distance=0.75),
689 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
690 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
691 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
694 just=c("left","top"),
695 gp=gpar(fontsize=48),
696 y=unit(1,"npc")-unit(0.25,"lines"),
697 x=unit(0,"npc")+unit(0.25,"lines"))
700 \caption{Change in charge forward ($ch_\mathrm{f}$) (A) or activation energy in
701 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the charge of the
702 monomer entering ($ch_\mathrm{monomer}$) and the average charge of
703 the vesicle ($\left<ch_\mathrm{vesicle}\right>$).}
704 \label{fig:chf_graph}
707 \subsubsection{Curvature Forward}
709 Curvature is a measure of the intrinsic propensity of specific lipids
710 to form micelles (positive curvature), inverted micelles (negative
711 curvature), or planar sheets (neutral
712 curvature)~\citep{Israelachvili1975:amphiphile_self_assembly}. %
713 In this formalism, curvature is measured as the ratio of the volume of
714 the lipid to the area of the head times the length of the lipid
715 ($S=\frac{v}{l_c a}$), so negative curvature is bounded by $(0,1)$,
716 neutral curvature is 1, and positive curvature is bounded by
717 $(1,\infty)$. The curvature can be transformed using $\log$, which has
718 the property of making the range of positive and negative curvature
719 equal, and distributed about 0.
721 As in the case of unsaturation, void formation is increased by the
722 presence of lipids with mismatched curvature. Thus, a larger
723 distribution of curvature in the vesicle increases the rate of lipid
724 insertion into the vesicle. However, a component with curvature
725 $cu^{-1}$ will cancel out a component with curvature $cu$, so we have
726 to log transform (turning these into $-\log cu$ and $\log cu$), then
727 take the absolute value ($\log cu$ and $\log cu$), and finally measure
728 the width of the distribution, which in the case of exactly mismatched
729 curvatures is $0$. Thus, by using the log transform to make the range
730 of the lipid curvature equal between positive and negative, and taking
731 the average to cancel out exactly mismatched curvatures, we come to an
732 equation with the shape $a^{\left<\log cu_\mathrm{vesicle}\right>}$.
733 An ad hoc base for $a$ which for the common case is $10$, yielding:
737 % cu_\mathrm{f} = 10^{\mathrm{stdev}\left|\log cu_\mathrm{vesicle}\right|}
738 cu_\mathrm{f} = 10^{\left|\left<\log cu_\mathrm{monomer} \right>\right|\mathrm{stdev} \left|\log cu_\mathrm{vesicle}\right|}
739 \label{eq:curvature_forward}
742 The most common $\left|\left<\log {cu}_v\right>\right|$ is around
743 $0.013$, which with the most common $\mathrm{stdev} \log
744 cu_\mathrm{vesicle}$ of $0.213$ leads to a $\Delta \Delta G^\ddagger$
745 of $\Sexpr{format(digits=3,to.kcal(10^(0.13*0.213)))}
746 \frac{\mathrm{kcal}}{\mathrm{mol}}$. This is a consequence of the
747 relatively matched curvatures in our environment. The full range of
748 $cu_\mathrm{f}$ values possible are shown in Figure~\ref{fig:cuf_graph}.
750 % 1.5 to 0.75 3 to 0.33
752 <<cuf_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
753 trellis.device(new=F,color=TRUE)
754 trellis.par.set(list(axis.line =list(col="transparent")))
755 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
756 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
757 grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
758 sd(abs(log(c(1,0.33)))),
759 sd(abs(log(c(0.33,3)))))),length.out=20),
760 y=seq(0,max(c(mean(log(c(1,3)),
761 mean(log(c(1,0.33))),
762 mean(log(c(0.33,3)))))),length.out=20))
763 grid$z <- 10^(grid$x*grid$y)
764 print(wireframe(z~x*y,grid,cuts=50,
767 scales=list(arrows=F,col=1,distance=0.75),
768 xlab=list(label=expression(paste(stdev,group("(",italic(log(cu[vesicle])),")"))),rot=30),
769 ylab=list(label=expression(paste(symbol("\341"),italic(log(cu[vesicle])),symbol("\361"))),rot=-35),
770 zlab=list(label=expression(italic(cu[f])),rot=93)),newpage=FALSE)
773 just=c("left","top"),
774 gp=gpar(fontsize=48),
775 y=unit(1,"npc")-unit(0.25,"lines"),
776 x=unit(0,"npc")+unit(0.25,"lines"))
778 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
779 grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
780 sd(abs(log(c(1,0.33)))),sd(abs(log(c(0.33,3)))))),length.out=20),
781 y=seq(0,max(c(mean(log(c(1,3)),
782 mean(log(c(1,0.33))),
783 mean(log(c(0.33,3)))))),length.out=20))
784 grid$z <- to.kcal(10^(grid$x*grid$y))
785 print(wireframe(z~x*y,grid,cuts=50,
788 scales=list(arrows=F,col=1,distance=0.75),
789 xlab=list(label=expression(paste(stdev,group("(",italic(log(cu[vesicle])),")"))),rot=30),
790 ylab=list(label=expression(paste(symbol("\341"),italic(log(cu[vesicle])),symbol("\361"))),rot=-35),
791 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
794 just=c("left","top"),
795 gp=gpar(fontsize=48),
796 y=unit(1,"npc")-unit(0.25,"lines"),
797 x=unit(0,"npc")+unit(0.25,"lines"))
800 \caption{Change in curvature forward ($cu_\mathrm{f}$) (A) or activation energy
801 in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the standard
802 deviation of the log of the curvature of the vesicle
803 ($\mathrm{stdev}\left(\log cu_\mathrm{vesicle}\right)$) and the mean
804 of the log of the curvature of the vesicle ($\left<\log
805 cu_\mathrm{vesicle}\right>$).}
806 \label{fig:cuf_graph}
809 \subsubsection{Length Forward}
811 As in the case of unsaturation, void formation is easier when vesicles
812 are made up of components of mismatched lengths. Thus, when the width
813 of the distribution of lengths is larger, the forward rate should be
814 greater as well, leading us to an equation of the form
815 $x^{\mathrm{stdev} l_\mathrm{vesicle}}$, where $\mathrm{stdev}
816 l_\mathrm{vesicle}$ is the standard deviation of the length of the
817 components of the vesicle, which has a maximum possible value of
818 $\Sexpr{format(digits=3,sd(c(rep(12,50),rep(24,50))))}$ and a minimum
819 of 0 in this set of simulations. Based on activation energy, a
820 reasonable base for $x$ is 2, leading to:
823 l_\mathrm{f} = 2^{\mathrm{stdev} l_\mathrm{vesicle}}
824 \label{eq:length_forward}
827 The most common $\mathrm{stdev} l_\mathrm{vesicle}$ is around $3.4$, which leads to
828 a range of $\Delta \Delta G^\ddagger$ of
829 $\Sexpr{format(digits=3,to.kcal(2^(3.4)))}
830 \frac{\mathrm{kcal}}{\mathrm{mol}}$.
832 % While it could be argued that increased length of the monomer could
833 % affect the rate of insertion into the membrane, it's not clear whether
834 % it would increase (by decreasing the number of available hydrogen
835 % bonds, for example) or decrease (by increasing the time taken to fully
836 % insert the acyl chain, for example) the rate of insertion or to what
837 % degree, so we do not take it into account in this formalism.
839 % \DLA{Incorporate McLean84 here}
840 % From McLean84LIB: Although it is difficult to measure cmc values for
841 % the sparingly soluble lipids used in this study, estimates for
842 % lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds,
843 % 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution;
844 % Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6
845 % X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X
846 % 10-8 M for DMPC was estimated from the linear relationship between ln
847 % cmc and the number of carbons in the PC acyl chain by using data for n
848 % = 7, 8, 10, and 16 [summarized in Tanford (1980)].
850 % From Nichols85: The magnitude of the dissociation rate constant
851 % decreases by a factor of approximately 3.2 per carbon increase in acyl
852 % chain length (see Table II here) {take into account for the formula;
855 % From Nichols85: The magnitude of the partition coefficient increases
856 % with acyl chain length [Keq(M-C6-NBD-PC) = (9.8±2.1) X 106 M-1 and Keq
857 % (P-C6-NBD-PC) = (9.4±1.3) x 107 M-1. {take into account for the
858 % formula; rz 8/17/2010}.
860 % From Nichols85: The association rate constant is independent of acyl
861 % chain length. {take into account for the formula; rz 8/17/2010}.
866 <<lf_graph,echo=FALSE,results='hide',fig.width=14,fig.height=5>>=
867 layout(matrix(1:2,nrow=1,ncol=2))
868 curve(2^x,from=0,to=sd(c(12,24)),
869 xlab=expression(paste(stdev,group("(",italic(l[vesicle]),")"))),
870 ylab=expression(italic(l[f])),
873 vps <- baseViewports()
874 pushViewport(vps$figure)
876 just=c("left","top"),
877 gp=gpar(fontsize=48),
878 y=unit(1,"npc")-unit(0.25,"lines"),
879 x=unit(0,"npc")+unit(0.25,"lines"))
881 curve(to.kcal(2^x),from=0,to=sd(c(12,24)),
882 xlab=expression(paste(stdev,group("(",italic(l[vesicle]),")"))),
883 ylab="Activation Energy (kcal/mol)",
885 vps <- baseViewports()
886 pushViewport(vps$figure)
888 just=c("left","top"),
889 gp=gpar(fontsize=48),
890 y=unit(1,"npc")-unit(0.25,"lines"),
891 x=unit(0,"npc")+unit(0.25,"lines"))
894 \caption{Change in length forward ($l_\mathrm{f}$) (A) or activation energy in
895 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the standard
896 deviation of the length of the vesicle
897 ($\mathrm{stdev}\left(l_\mathrm{vesicle}\right)$).}
901 \subsubsection{Complex Formation}
902 There is no contribution of complex formation to the forward reaction
903 rate in the current formalism.
907 \label{eq:complex_formation_forward}
910 \subsection{Backward adjustments ($k_{\mathrm{b}i\mathrm{adj}}$)}
912 Just as the forward rate constant adjustment $k_{\mathrm{f}i\mathrm{adj}}$
913 does, the backwards rate constant adjustment $k_{\mathrm{b}i\mathrm{adj}}$
914 takes into account unsaturation ($un_\mathrm{b}$), charge ($ch_\mathrm{b}$), curvature
915 ($cu_\mathrm{b}$), length ($l_\mathrm{b}$), and complex formation ($CF1_\mathrm{b}$), each of
916 which are modified depending on the specific component and the vesicle
917 from which the component is exiting:
921 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}
925 \subsubsection{Unsaturation Backward}
927 Unsaturation also influences the ability of a lipid molecule to leave
928 a membrane. If a molecule has an unsaturation level which is different
929 from the surrounding membrane, it will be more likely to leave the
930 membrane. The more different the unsaturation level is, the greater
931 the propensity for the lipid molecule to leave. However, a vesicle
932 with some unsaturation (eg. average unsaturation of 2) is more
933 favorable for lipids with more unsaturation (eg. 3) than lipids with
934 equivalently less unsaturation (eg. 1), so the difference in energy
935 between unsaturation is not linear. Therefore, an equation with the
936 shape $x^{\left| y^{-\left< un_\mathrm{vesicle}\right>
937 }-y^{-un_\mathrm{monomer}} \right| }$, where
938 $\left<un_\mathrm{vesicle}\right>$ is the average unsaturation of the
939 vesicle and $un_\mathrm{monomer}$ is the unsaturation of the monomer,
940 allows for increasing the efflux of molecules from membranes where
941 they strongly mismatch, while allowing vesicles with greater
942 unsaturation to tolerate greater unsaturation mismatch between
943 monomers and the vesicle. The bases x and y were chosen ad hoc to
944 produce reasonable Eyring energies of activation over the range of
945 unsaturations expected, leading to:
949 un_\mathrm{b} = 7^{1-\left(20\left(2^{-\left<un_\mathrm{vesicle} \right>} - 2^{-un_\mathrm{monomer}} \right)^2+1\right)^{-1}}
950 \label{eq:unsaturation_backward}
953 The most common $\left<un_\mathrm{vesicle}\right>$ is around $1.7$, which
954 leads to a range of $\Delta \Delta G^\ddagger$ from
955 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-0)^2+1))))}
956 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with 0 unsaturation
958 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-4)^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
959 for monomers with 4 unsaturations. See Figure~\ref{fig:unb_graph} for
960 the full range of possible values.
964 <<unb_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
965 trellis.device(new=F,color=TRUE)
966 trellis.par.set(list(axis.line =list(col="transparent")))
967 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
968 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
969 grid <- expand.grid(x=seq(0,4,length.out=20),
970 y=seq(0,4,length.out=20))
971 grid$z <- (7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1)))
972 print(wireframe(z~x*y,grid,cuts=50,
975 scales=list(arrows=F,col=1,distance=0.75),
976 xlab=list(label=expression(paste(symbol("\341"),italic(un[vesicle]),symbol("\361"))),rot=30),
977 ylab=list(label=expression(italic(un[monomer])),rot=-35),
978 zlab=list(label=expression(italic(un[b])),rot=93)),newpage=FALSE)
980 just=c("left","top"),
981 gp=gpar(fontsize=48),
982 y=unit(1,"npc")-unit(0.25,"lines"),
983 x=unit(0,"npc")+unit(0.25,"lines"))
985 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
986 grid <- expand.grid(x=seq(0,4,length.out=20),
987 y=seq(0,4,length.out=20))
988 grid$z <- to.kcal((7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1))))
989 print(wireframe(z~x*y,grid,cuts=50,
992 scales=list(arrows=F,col=1,distance=0.75),
993 xlab=list(label=expression(paste(symbol("\341"),italic(un[vesicle]),symbol("\361"))),rot=30),
994 ylab=list(label=expression(italic(un[monomer])),rot=-35),
995 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
998 just=c("left","top"),
999 gp=gpar(fontsize=48),
1000 y=unit(1,"npc")-unit(0.25,"lines"),
1001 x=unit(0,"npc")+unit(0.25,"lines"))
1004 \caption{Change in unsaturation backward ($un_\mathrm{b}$) (A) or activation energy in
1005 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the unsaturation of the
1006 monomer leaving ($un_\mathrm{monomer}$) and the average unsaturation of
1007 the vesicle ($\left<un_\mathrm{vesicle}\right>$).}
1008 \label{fig:unb_graph}
1011 \subsubsection{Charge Backwards}
1012 As in the case of monomers entering a vesicle, monomers leaving a
1013 vesicle leave faster if their charge has the same sign as the average
1014 charge vesicle. An equation of the form $ch_\mathrm{b} = a^{\left<ch_v\right>
1015 ch_m}$ is then appropriate, and using a base of $a=20$ yields:
1018 ch_\mathrm{b} = 20^{\left<{ch}_v\right> {ch}_m}
1019 \label{eq:charge_backwards}
1022 The most common $\left<ch_v\right>$ is around $-0.164$, which leads to
1023 a range of $\Delta \Delta G^\ddagger$ from
1024 $\Sexpr{format(digits=3,to.kcal(20^(-.164*-1)))}
1025 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with charge $-1$ to
1026 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with charge $0$.
1027 See Figure~\ref{fig:chb_graph} for the full range of possible values
1032 <<chb_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1033 trellis.device(new=F,color=TRUE)
1034 trellis.par.set(list(axis.line =list(col="transparent")))
1035 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1036 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1037 x <- seq(-1,0,length.out=20)
1038 y <- seq(-1,0,length.out=20)
1039 grid <- expand.grid(x=x,y=y)
1040 grid$z <- as.vector(20^(outer(x,y)))
1041 print(wireframe(z~x*y,grid,cuts=50,
1044 scales=list(arrows=F,col=1,distance=0.75),
1045 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
1046 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
1047 zlab=list(label=expression(italic(ch[b])),rot=93)),newpage=FALSE)
1050 just=c("left","top"),
1051 gp=gpar(fontsize=48),
1052 y=unit(1,"npc")-unit(0.25,"lines"),
1053 x=unit(0,"npc")+unit(0.25,"lines"))
1055 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1056 x <- seq(-1,0,length.out=20)
1057 y <- seq(-1,0,length.out=20)
1058 grid <- expand.grid(x=x,y=y)
1059 grid$z <- to.kcal(as.vector(20^(outer(x,y))))
1060 print(wireframe(z~x*y,grid,cuts=50,
1063 scales=list(arrows=F,col=1,distance=0.75),
1064 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
1065 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
1066 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1069 just=c("left","top"),
1070 gp=gpar(fontsize=48),
1071 y=unit(1,"npc")-unit(0.25,"lines"),
1072 x=unit(0,"npc")+unit(0.25,"lines"))
1075 \caption{Change in charge backward ($ch_\mathrm{b}$) (A) or activation energy in
1076 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the charge of the
1077 monomer leaving ($ch_\mathrm{monomer}$) and the average charge of
1078 the vesicle ($\left<ch_\mathrm{vesicle}\right>$).}
1079 \label{fig:chb_graph}
1083 \subsubsection{Curvature Backwards}
1085 The less a monomer's intrinsic curvature matches the average curvature
1086 of the vesicle in which it is in, the greater its rate of efflux. If
1087 the difference is 0, $cu_\mathrm{f}$ needs to be one. To map negative and
1088 positive curvature to the same range, we also need take the logarithm.
1089 Positive and negative curvature lipids cancel each other out,
1090 resulting in an average curvature of 0; the average of the log also
1091 has this property. Increasing mismatches in curvature increase the
1092 rate of efflux, but asymptotically. An equation which satisfies these
1093 criteria has the form $cu_\mathrm{f} = a^{1-\left(b\left( \left< \log
1094 cu_\mathrm{vesicle} \right> -\log
1095 cu_\mathrm{monomer}\right)^2+1\right)^{-1}}$. An alternative
1096 form would use the absolute value of the difference, however, this
1097 yields a cusp and sharp increase about the curvature equilibrium. We
1098 have chosen bases of $a=7$ and $b=20$ ad hoc, based on activation
1099 energy considerations.
1102 cu_\mathrm{b} = 7^{1-\left(20\left(\left< \log cu_\mathrm{vesicle} \right> -\log cu_\mathrm{monomer} \right)^2+1\right)^{-1}}
1103 \label{eq:curvature_backwards}
1106 The most common $\left<\log cu_\mathrm{vesicle}\right>$ is around
1107 $-0.013$, which leads to a range of $\Delta \Delta G^\ddagger$ from
1108 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(0.8))^2+1))))}
1109 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature 0.8 to
1110 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(1.3))^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
1111 for monomers with curvature 1.3 to
1112 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature near
1113 1. The full range of values possible for $cu_\mathrm{b}$ are shown in
1114 Figure~\ref{fig:cub_graph}
1116 % \RZ{What about the opposite curvatures that actually do fit to each
1120 % From Kumar91: It has also been shown experimentally that molecules
1121 % having complementary molecular shapes can form bilayer structures.
1122 % Inverted cone (wedge)- shaped lysoPC and cone-shaped Chol or
1123 % unsaturated PE can interact stoichiometrically to form bilayer
1124 % structures (9-12). {quote to justify our formula for complementary
1129 <<cub_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1130 trellis.device(new=F,color=TRUE)
1131 trellis.par.set(list(axis.line =list(col="transparent")))
1132 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1133 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1134 grid <- expand.grid(x=log(seq(0.8,1.33,length.out=20)),
1135 y=seq(0.8,1.33,length.out=20))
1136 grid$z <- 7^(1-1/(20*(grid$x-log(grid$y))^2+1))
1137 print(wireframe(z~x*y,grid,cuts=50,
1140 scales=list(arrows=F,col=1,distance=0.75),
1141 xlab=list(label=expression(paste(symbol("\341"),log(italic(cu[vesicle])),symbol("\361"))),rot=30),
1142 ylab=list(label=expression(italic(cu[monomer])),rot=-35),
1143 zlab=list(label=expression(italic(cu[b])),rot=93)),newpage=FALSE)
1146 just=c("left","top"),
1147 gp=gpar(fontsize=48),
1148 y=unit(1,"npc")-unit(0.25,"lines"),
1149 x=unit(0,"npc")+unit(0.25,"lines"))
1151 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1152 grid <- expand.grid(x=log(seq(0.8,1.33,length.out=20)),
1153 y=seq(0.8,1.33,length.out=20))
1154 grid$z <- to.kcal(7^(1-1/(20*(grid$x-log(grid$y))^2+1)))
1155 print(wireframe(z~x*y,grid,cuts=50,
1158 scales=list(arrows=F,col=1,distance=0.75),
1159 xlab=list(label=expression(paste(symbol("\341"),log(italic(cu[vesicle])),symbol("\361"))),rot=30),
1160 ylab=list(label=expression(italic(cu[monomer])),rot=-35),
1161 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1164 just=c("left","top"),
1165 gp=gpar(fontsize=48),
1166 y=unit(1,"npc")-unit(0.25,"lines"),
1167 x=unit(0,"npc")+unit(0.25,"lines"))
1170 \caption{Change in curvature backward ($cu_\mathrm{b}$) (A) or activation
1171 energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the
1172 curvature of the monomer leaving ($ch_\mathrm{monomer}$) and the
1173 average of the log of the curvature of the vesicle
1174 ($\left<\log cu_\mathrm{vesicle}\right>$).}
1175 \label{fig:cub_graph}
1179 \subsubsection{Length Backwards}
1181 In a model membrane, the dissociation constant increases by a factor
1182 of approximately 3.2 per carbon decrease in acyl chain length~%
1183 \citep{Nichols1985:phospholipid_monomer_vesicle_thermodynamics}.
1184 Unfortunately, the experimental data known to us only measures the
1185 effect of chain lengths less than or equal to the bulk lipid, not the
1186 effect of lengths exceeding it, and is only known for one bulk lipid
1187 component (DOPC). We assume then, that the increase is in relationship
1188 to the average vesicle, and that lipids with larger acyl chain length
1189 will also show an increase in their dissociation constant.
1192 l_\mathrm{b} = 3.2^{\left|\left<l_\mathrm{vesicle}\right>-l_\mathrm{monomer}\right|}
1193 \label{eq:length_backward}
1196 The most common $\left<l_\mathrm{vesicle}\right>$ is around $17.75$,
1197 which leads to a range of $\Delta \Delta G^\ddagger$ from
1198 $\Sexpr{format(digits=3,to.kcal(3.2^abs(12-17.75)))}
1199 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with length 12 to
1200 $\Sexpr{format(digits=3,to.kcal(3.2^abs(24-17.75)))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
1201 for monomers with length 24 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$
1202 for monomers with length near 18. The full range of possible values of
1203 $l_\mathrm{b}$ are shown in Figure~\ref{fig:lb_graph}
1205 % (for methods? From McLean84LIB: The activation free energies and free
1206 % energies of transfer from self-micelles to water increase by 2.2 and
1207 % 2.1 kJ mol-' per methylene group, respectively. \RZ{see if we can use it
1208 % to justify arranging our changed in activating energy around 1
1213 <<lb_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1214 trellis.device(new=F,color=TRUE)
1215 trellis.par.set(list(axis.line =list(col="transparent")))
1216 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1217 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1218 grid <- expand.grid(x=seq(12,24,length.out=20),
1219 y=seq(12,24,length.out=20))
1220 grid$z <- 3.2^(abs(grid$x-grid$y))/1e6
1221 print(wireframe(z~x*y,grid,cuts=50,
1224 scales=list(arrows=F,col=1,distance=0.75),
1225 xlab=list(label=expression(paste(symbol("\341"),italic(l[vesicle]),symbol("\361"))),rot=30),
1226 ylab=list(label=expression(italic(l[monomer])),rot=-35),
1227 zlab=list(label=expression(italic(l[b])%*%10^-6),rot=93)),newpage=FALSE)
1230 just=c("left","top"),
1231 gp=gpar(fontsize=48),
1232 y=unit(1,"npc")-unit(0.25,"lines"),
1233 x=unit(0,"npc")+unit(0.25,"lines"))
1235 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1236 grid <- expand.grid(x=seq(12,24,length.out=20),
1237 y=seq(12,24,length.out=20))
1238 grid$z <- to.kcal(3.2^(abs(grid$x-grid$y)))
1239 print(wireframe(z~x*y,grid,cuts=50,
1242 scales=list(arrows=F,col=1,distance=0.75),
1243 xlab=list(label=expression(paste(symbol("\341"),italic(l[vesicle]),symbol("\361"))),rot=30),
1244 ylab=list(label=expression(italic(l[monomer])),rot=-35),
1245 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1248 just=c("left","top"),
1249 gp=gpar(fontsize=48),
1250 y=unit(1,"npc")-unit(0.25,"lines"),
1251 x=unit(0,"npc")+unit(0.25,"lines"))
1254 \caption{Change in length backwards ($l_\mathrm{b}\times 10^{-6}$) (A) or activation energy
1255 in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the length of the
1256 monomer leaving ($l_\mathrm{monomer}$) and the average length of
1257 lipids in the vesicle ($\left<l_\mathrm{vesicle}\right>$).}
1258 \label{fig:lb_graph}
1262 \subsubsection{Complex Formation Backward}
1264 Complex formation ($CF1$) describes the interaction between CHOL and
1265 PC or SM, where PC or SM protects the hydroxyl group of CHOL from
1266 interactions with water %
1267 \citep{Huang1999:chol_solubility_model,
1268 Huang1999:cholesterol_solubility, McConnell2006, McConnell2003}%
1269 . PC ($CF1=2$) can interact with two CHOL, and SM ($CF1=3$) with three
1270 CHOL ($CF1=-1$). If the average of $CF1$ is positive (excess of SM and
1271 PC with regards to complex formation), components with negative $CF1$
1272 (CHOL) will be retained. If average $CF1$ is negative, components with
1273 positive $CF1$ are retained. An equation which has this property is
1274 $CF1_\mathrm{b}=a^{\left<CF1_\mathrm{vesicle}\right>
1275 CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{vesicle}\right>
1276 CF1_\mathrm{monomer}\right|}$, where difference of the exponent is
1277 zero if the average $CF1$ and the $CF1$ of the component have the same
1278 sign, or double the product if the signs are different. Based on
1279 activation energy considerations, we took an ad hoc base for $a$ as
1284 CF1_\mathrm{b}=1.5^{\left<CF1_\mathrm{vesicle}\right> CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{vesicle}\right> CF1_\mathrm{monomer}\right|}
1285 \label{eq:complex_formation_backward}
1288 The most common $\left<CF1_\mathrm{vesicle}\right>$ is around $0.925$,
1289 which leads to a range of $\Delta \Delta G^\ddagger$ from
1290 $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*-1-abs(0.925*-1))))}
1291 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
1293 $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*2-abs(0.925*2))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
1294 for monomers with complex formation $2$ to
1295 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
1296 formation $0$. The full range of possible values for $CF1_\mathrm{b}$ are
1297 depicted in Figure~\ref{fig:cf1b_graph}.
1302 <<cf1b_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1303 trellis.device(new=F,color=TRUE)
1304 trellis.par.set(list(axis.line =list(col="transparent")))
1305 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1306 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1307 grid <- expand.grid(x=seq(-1,3,length.out=20),
1308 y=seq(-1,3,length.out=20))
1309 grid$z <- 1.5^(grid$x*grid$y-abs(grid$x*grid$y))
1310 print(wireframe(z~x*y,grid,cuts=50,
1313 scales=list(arrows=F,col=1,distance=0.75),
1314 xlab=list(label=expression(paste(symbol("\341"),italic(CF1[vesicle]),symbol("\361"))),rot=30),
1315 ylab=list(label=expression(italic(CF1[monomer])),rot=-35),
1316 zlab=list(label=expression(italic(CF1[b])),rot=93)),newpage=FALSE)
1319 just=c("left","top"),
1320 gp=gpar(fontsize=48),
1321 y=unit(1,"npc")-unit(0.25,"lines"),
1322 x=unit(0,"npc")+unit(0.25,"lines"))
1324 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1325 grid <- expand.grid(x=seq(-1,3,length.out=20),
1326 y=seq(-1,3,length.out=20))
1327 grid$z <- to.kcal(1.5^(grid$x*grid$y-abs(grid$x*grid$y)))
1328 print(wireframe(z~x*y,grid,cuts=50,
1331 scales=list(arrows=F,col=1,distance=0.75),
1332 xlab=list(label=expression(paste(symbol("\341"),italic(CF1[vesicle]),symbol("\361"))),rot=30),
1333 ylab=list(label=expression(italic(CF1[monomer])),rot=-35),
1334 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1337 just=c("left","top"),
1338 gp=gpar(fontsize=48),
1339 y=unit(1,"npc")-unit(0.25,"lines"),
1340 x=unit(0,"npc")+unit(0.25,"lines"))
1343 \caption{Change in complex formation 1 backwards ($CF1_\mathrm{b}$) (A) or
1344 activation energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus
1345 the complex formation of the monomer leaving
1346 ($CF1_\mathrm{monomer}$) and the average complex formation of the
1347 vesicle ($\left<CF1_\mathrm{vesicle}\right>$).}
1348 \label{fig:cf1b_graph}
1351 \section{Simulation Methodology}
1353 \subsection{Overall Architecture}
1355 The simulations are currently run using a program written in perl with
1356 various modules to handle the subsidiary parts. It produces output for
1357 each generation, including the step immediately preceding and
1358 immediately following a vesicle split (and optionally, each step) that
1359 is written to a state file which contains the state of the vesicle,
1360 environment, kinetic parameters, program invocation options, time, and
1361 various other parameters necessary to recreate the state vector at a
1362 given time. This output file is then read by a separate program in
1363 perl to produce different output which is generated out-of-band;
1364 output which includes graphs and statistical analysis is performed
1365 using R~\citep{R:maincite} (and various grid graphics modules) called
1366 from the perl program.
1368 The separation of simulation and output generation allows refining
1369 output, and simulations performed on different versions of the
1370 underlying code to be compared using the same analysis methods and
1371 code. It also allows later simulations to be restarted from a specific
1372 generation, utilizing the same environment.
1374 Random variables of different distributions are calculated using
1375 Math::Random (0.71), which is seeded for each run using entropy from
1376 the linux kernel's urandom device.
1378 Code is available upon request.
1380 \subsection{Environment Creation}
1382 \subsubsection{Components}
1383 The environment contains concentrations of components. In the current
1384 set of simulations, there are \Sexpr{1+4*5*7} different possible
1385 components, consisting of PC, PE, PS, SM, and CHOL, with all lipids
1386 except for CHOL having 5 possible unsaturations ranging from 0 to 4,
1387 and 7 lengths from $12,14,...,22$ ($7\cdot
1388 5\cdot4+1=\Sexpr{1+4*5*7}$). In cases where the environment has less
1389 than the maximum number of components, the components are selected in
1390 order without replacement from a randomly shuffled deck of components
1391 (with the exception of CHOL) represented once until the desired number
1392 of unique components is obtained. CHOL is overrepresented
1393 $\Sexpr{5*7}$ times to be at the level of other lipid types, ensuring
1394 that the probability of CHOL being absent in the environment is the
1395 same as the probability of any one of the other lipid types (PS, PE,
1396 etc.) being absent. This reduces the number of simulations with a
1397 small number of components which are devoid of CHOL.
1399 \subsubsection{Concentrations}
1400 Once the components of the environment have been selected, their
1401 concentrations are determined. In experiments where the environmental
1402 concentration is the same across all lipid components, the
1403 concentration is set to $10^{-10}\mathrm{M}$. In other cases, the
1404 environmental concentration is set to a random number from a gamma
1405 distribution with shape parameter 2 and an average of
1406 $10^{-10}\mathrm{M}$. The base concentration ($10^{-10}\mathrm{M}$)
1407 can also be altered at the initialization of the experiment to
1408 specific values for specific lipid types or components.
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. Lipid molecules are then selected to be
1417 added to the lipid membrane according to four different methods. In
1418 the constant method, molecules are added in direct proportion to their
1419 concentration in the environment. The uniform method adds molecules in
1420 proportion to their concentration in the environment scaled by a
1421 uniform random value, whereas the random method adds molecules in
1422 proportion to a uniform random value. The final method is a binomial
1423 method, which adjust the probability of adding a molecule of a
1424 specific component by the concentration of that component, and then
1425 adds the molecules one by one to the membrane. This final method is
1426 also used in the first three methods to add any missing molecules to
1427 the starting vesicle which were unallocated due to the requirement to
1428 add integer numbers of molecules. (For example, if a vesicle was to
1429 have 10 molecules split evenly between three lipid types, the 10th
1430 molecule would be assigned by randomly choosing between the three
1431 lipid types, weighted by their concentration in the environment.)
1433 \subsection{Simulation Step}
1435 Once the environment has been created and the initial vesicle has been
1436 formed, molecules join and leave the vesicle based on the kinetic
1437 parameters and state equation discussed above until the vesicle splits
1438 forming two progeny vesicles, one of which the program continues to
1441 \subsubsection{Calculation of Vesicle Properties}
1443 \label{sec:ves_prop}
1444 $S_\mathrm{vesicle}$ is the surface area of the vesicle, and is the sum of
1445 the surface area of all of the individual lipid components; each lipid
1446 type has a different surface area; we assume that the lipid packing is
1447 optimal, and there is no empty space.
1449 For the other vesicle properties (length, unsaturation, charge, and
1450 curvature), we calculate the mean, the standard deviation, the mean of
1451 the log, the mean of the absolute value of the log, the standard
1452 deviation of the log, and the standard deviation of the absolute value
1453 of the log. All cases, when we take the log, we take the log of the
1454 absolute value, and map $\log(0)$ to $0$. For the purpose of plotting
1455 vesicle properties, we only plot the mean of the property.
1457 \subsubsection{Joining and Leaving of Lipid Molecules}
1459 Determining the number of molecules to add to the lipid membrane
1460 ($n_i$) requires knowing $k_{\mathrm{f}i_\mathrm{adj}}$, the surface area of the
1461 vesicle $S_\mathrm{vesicle}$ (see Section \ref{sec:ves_prop}), the time interval
1462 $dt$ during which lipids are added, the base $k_{\mathrm{f}i}$, and the
1463 concentration of the monomer in the environment
1464 $\left[C_{i\mathrm{vesicle}}\right]$ (see Equation~\ref{eq:state}).
1465 $k_{\mathrm{f}i\mathrm{adj}}$ is calculated (see Equation~\ref{eq:kf_adj}) based on the
1466 vesicle properties and their hypothesized effect on the rate (in as
1467 many cases as possible, experimentally based)
1468 (see Section~\ref{sec:kinetic_adjustments} for details). $dt$ can be varied
1469 (see Section~\ref{sec:step_duration}), but for a given step is constant. This
1470 leads to the following:
1472 $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$
1474 In the cases where $n_i > 1$, the integer number of molecules is
1475 added. Fractional $n_i$ or the fractional remainder after the addition
1476 of the integer molecules are the probability of adding a specific
1477 molecule, and are compared to a uniformly distributed random value
1478 between 0 and 1. If the random value is less than or equal to the
1479 fractional part of $n_i$, an additional molecule is added.
1481 Molecules leaving the vesicle are handled in a similar manner, with
1483 $n_i = k_{\mathrm{b}i}k_{\mathrm{b}i_\mathrm{adj}}C_{i_\mathrm{vesicle}}\mathrm{N_A}dt$.
1485 Before addition or removal of molecules, the properties of the vesicle
1486 are calculated; this is done so that molecules can be considered to be
1487 added and removed at the same instant, even though additions are
1488 handled first programmatically. This also avoids cases where a removal
1489 would have resulted in a negative number of molecules for a particular
1492 \subsubsection{Step duration}
1493 \label{sec:step_duration}
1495 $dt$ is the time taken for each step of joining, leaving, and checking
1496 split. In the current implementation, it starts with a value of
1497 $10^{-6}$~$\mathrm{s}$ but this is modified in between each step if the
1498 number of molecules joining or leaving is too large or small. If more
1499 than half of the starting vesicle size molecules join or leave in a
1500 single step, $dt$ is reduced by half. If less than the starting
1501 vesicle size molecules join or leave in 100 steps, $dt$ is doubled.
1502 This is necessary to curtail run times and to automatically adjust to
1503 different experimental runs.
1505 \subsubsection{Vesicle splitting}
1507 If a vesicle has grown to a size which is more than double the
1508 starting vesicle size, the vesicle splits. Molecules are assigned to
1509 the progeny vesicles at random, with each progeny vesicle having an
1510 equal chance of getting a single molecule. The number of molecules to
1511 assign to the first vesicle has binomial distribution with a
1512 probability of an event in each trial of 0.5 and a number of trials
1513 equal to the number of molecules.
1517 The environment, initial vesicle, and the state of the vesicle
1518 immediately before and immediately after splitting are stored
1519 to produce later output.
1521 \section{Analyzing output}
1523 The analysis of output is handled by a separate perl program which
1524 shares many common modules with the simulation program. Current output
1525 includes simulation progress, summary tables, summary statistics, and
1528 \subsection{PCA plots}
1530 Two major groups of axes that contribute to vesicle variation between
1531 subsequent generations are the components and properties of vesicles.
1532 Each component in a vesicle is an axis on its own; it can be measured
1533 either as an absolute number of molecules in each component, or the
1534 fraction of molecules of that component over the total number of
1535 molecules; the second approach is often more convenient, as it allows
1536 vesicles with different number of molecules to be directly compared
1537 (though it hides any effect of vesicle size).
1539 In order to visualize the variation between and transition of
1540 subsequent generations of vesicles from their initial state in the
1541 simulation, to their final state at the termination of the simulation,
1542 we plot the projection of the generations onto the two principle PCA
1543 axes (and alternatively, any pairing of the axes).
1545 \subsection{Carpet plots}
1547 Carpet plots show the distance/similarity between the composition or
1548 properties of all generations in a single run. The difference between
1549 each group of vesicle can be calculated using Euclidean distance
1550 (properties and compositions) or H similarity (composition only). We
1551 must use Euclidean distance for properties because the H distance
1552 metric is invalid when comparing negative vector elements to positive
1555 In addition to showing distance or similarity, carpet plots also
1556 depict propersomes and composomes as square boxes on the diagonals and
1557 propertypes and compotypes as rectangles off the diagonals, each
1558 propertype or compotype with a distinct color.
1560 \subsection{Propersomes, propertypes, composomes and compotypes}
1562 A generation is considered to be a propersome if it is less than $0.6$
1563 units (by Euclidean distance of normalized properties) away from the
1564 generation immediately following and preceding. Likewise, a generation
1565 is in a composome if its H similarity is more than $0.9$ (by the
1566 normalized H metric) from the preceding and following generations.
1567 Propersomes and composomes are then assigned to propertypes and
1568 compotypes using Paritioning Around Medioids
1569 (PAM). All values of $k$ between 2 and 15
1570 (or the number of propersomes and composomes, if that is less than 15)
1571 are tried, and the clustering with the smallest
1572 silhouette~\citep{Rousseeuw1987:silhouettes} is chosen as the ideal
1573 clustering~\citep{Shenhav2005:pgard}.
1576 \bibliographystyle{unsrtnat}
1577 \bibliography{references.bib}