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[noblocks]{authblk}
26 \usepackage[backend=biber,natbib=true,hyperref=true,citestyle=numeric-comp,style=nature,autocite=inline]{biblatex}
27 \addbibresource{references.bib}
28 \usepackage[hyperfigures,bookmarks,colorlinks,citecolor=black,filecolor=black,linkcolor=black,urlcolor=black]{hyperref}
29 \usepackage[capitalise]{cleveref}
30 %\usepackage[sectionbib,sort&compress,numbers]{natbib}
31 % \usepackage[nomargin,inline,draft]{fixme}
32 %\usepackage[x11names,svgnames]{xcolor}
33 %\usepackage{texshade}
34 \renewcommand{\textfraction}{0.15}
35 \renewcommand{\topfraction}{0.85}
36 \renewcommand{\bottomfraction}{0.65}
37 \renewcommand{\floatpagefraction}{0.60}
38 %\renewcommand{\baselinestretch}{1.8}
40 \newcommand{\DLA}[1]{\textcolor{red}{\fxnote{DLA: #1}}}
42 \newcommand{\RZ}[1]{\textcolor{red}{\fxnote{RZ: #1}}}
44 \newcommand{\OM}[1]{\textcolor{red}{\fxnote{OM: #1}}}
55 \title{Supplemental Material}
57 % rubber: module bibtex
58 % rubber: module natbib
63 <<load.libraries,echo=FALSE,results="hide",warning=FALSE,message=FALSE,error=FALSE,cache=FALSE>>=
64 opts_chunk$set(dev="CairoPDF",out.width="\\columnwidth",out.height="0.7\\textheight",out.extra="keepaspectratio")
65 opts_chunk$set(cache=TRUE, autodep=TRUE)
66 options(scipen = -2, digits = 1)
71 to.latex <- function(x){
72 gsub("\\\\","\\\\\\\\",latexSN(x))
75 to.kcal <- function(k,temp=300) {
77 return(-gasconst*temp*log(k)/1000)
81 % \setcounter{figure}{0} \setcounter{table}{0}
83 \renewcommand{\thefigure}{S\@arabic\c@figure}
84 \renewcommand{\thetable}{S\@arabic\c@table}
87 % \section{Competition Implementation}
88 % \subsection{Implementation changes}
91 % \item settable maximum number of vesicles to track (default $10^4$)
92 % \item start with 1~L ($10^{-3}$~m$^3$) cube
93 % \item if at any point the number of vesicles exceeds the maximum
94 % number, chop the volume and environment molecule number into tenths,
95 % randomly select one tenth of the vesicles, and continue tracking.
96 % \item generations will be counted per vesicle, and each progeny
97 % vesicle will have a generation number one greater than the parental
99 % \item 100 generations can result in as many as $2^{100}$
100 % ($\Sexpr{2^100}$) vesicles or as few as
102 % \item Environment will use a specific number of each component instead
103 % of a constant concentration; as the number may be larger than
104 % \texttt{long long} ($2^{64}$), we use libgmp to handle an arbitrary
105 % precision number of components
108 % \subsection{Infrastructure changes}
110 % \item Rewrite core bits in C
111 % \item Use libgmp for handling large ints
112 % \item Use openmpi to split the calculations out over multiple
113 % machines/processors and allow deploying to larger
114 % clusters/supercomputers
119 \section{State Equation}
120 % double check this with the bits in the paper
122 For a system with monomers $(_\mathrm{monomer})$ and a vesicle
123 $(_\mathrm{vesicle})$, the change in concentration of the $i$th component of
124 a lipid vesicle per change in time ($d \left[C_{i_\mathrm{vesicle}}\right]/dt$)
125 can be described by a modification of the basic mass action law:
128 \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] -
129 k_{\mathrm{b}i}k_{\mathrm{b}i\mathrm{adj}}\left[C_{i_\mathrm{vesicle}}\right]
133 The base forward kinetic parameter for the $i$th component is $k_{\mathrm{f}i}$
134 and is dependent on the particular lipid type (\ac{PC}, \ac{PS}, \ac{SM}, etc.). The
135 forward adjustment parameter, $k_{\mathrm{f}i\mathrm{adj}}$, is based on the
136 properties of the vesicle and the specific component (type, length,
137 unsaturation, etc.) (see \cref{eq:kf_adj,sec:kinetic_adjustments}).
138 $\left[C_{i_\mathrm{monomer}}\right]$ is the molar concentration of
139 monomer of the $i$th component. $\left[S_\mathrm{vesicle}\right]$ is the surface
140 area of the vesicle per volume. The base backwards kinetic parameter
141 for the $i$th component is $k_{\mathrm{b}i}$ and its adjustment parameter
142 $k_{\mathrm{b}i\mathrm{adj}}$ (see \cref{eq:kb_adj,sec:kinetic_adjustments}).
143 $\left[C_{i_\mathrm{vesicle}}\right]$ is the molar concentration of
144 the $i$th component in the vesicle.
146 \subsection{Per-Lipid Kinetic Parameters}
148 <<kf_prime,echo=FALSE,results='hide'>>=
149 kf.prime <- c(3.7e6,3.7e6,5.1e7,3.7e6,2.3e6)
150 kf <- (as.numeric(kf.prime)*10^-3)/(63e-20*6.022e23)
153 Each of the 5 lipid types has different kinetic parameters; where
154 available, these were taken from literature.
158 \begin{tabular}{c c c c c c c c}
160 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({Å}^2\right)$ & Charge & $\mathrm{CF}1$ & Curvature \\
162 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 \\
163 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 \\
164 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 \\
165 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 \\
166 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 \\
169 \caption{Kinetic parameters and molecular properties of lipid types}
170 \label{tab:kinetic_parameters_lipid_types}
173 %%% \DLA{I think we may just reduce these three sections; area, $k_\mathrm{f}$
174 %%% and $k_\mathrm{b}$ to \cref{tab:kinetic_parameters_lipid_types} with
177 %%% \RZ{Yes, but we also have to have then as comments the numbers that
178 %%% are not supported by refs}
180 \subsubsection{$k_\mathrm{f}$ for lipid types}
181 $k_{\mathrm{f}_\mathrm{PC}}$ was measured by
182 \citet{Nichols1985:phospholipid_monomer_vesicle_thermodynamics} and
183 was found to be $3.7\times 10^6$~$\frac{1}{\mathrm{M}\mathrm{s}}$ by the
184 partitioning of \ac{NBDPC} between \ac{DOPC} vesicles and water. As
185 similar references are not available for \ac{SM} or \ac{PS}, we assume that they have
186 the same $k_\mathrm{f}$. For \ac{CHOL}, no direct measurement of $k_\mathrm{f}$ is available,
187 however, \citet{Estronca2007:dhe_kinetics} measured the transfer of
188 \ac{DHE} from \ac{BSA} to \acp{LUV}, and found a $k_\mathrm{f}$ of $5.1\times 10^7$~%
189 $\frac{1}{\mathrm{M} \mathrm{s}}$. We assume that this value is close
190 to that of \ac{CHOL}, and use it for $k_{\mathrm{f}_\mathrm{\ac{CHOL}}}$. In the case of
191 \ac{PE}, \citet{Abreu2004:kinetics_ld_lo} measured the association of
192 \ac{NBDDMPE} with \ac{POPC} \acp{LUV} found a value for $k_\mathrm{f}$ of $2.3 \times 10^{6}$~%
193 $\frac{1}{\mathrm{M} \mathrm{s}}$. These three authors used a slightly
194 different kinetic formalism ($\frac{d\left[A_v\right]}{dt} =
195 k'_\mathrm{f}[A_m][L_v] - k_\mathrm{b}[A_v]$), so we correct their values of $k_\mathrm{f}$ by
196 multiplying by the surface area of a mole of lipids.
198 \subsubsection{$k_\mathrm{b}$ for lipid types}
200 \citet{Wimley1990:dmpc_exchange} measured the half time for the
201 exchange of [$^3$H]\ac{DMPC} between \acp{LUV} at 30\textdegree C, and found it
202 to be 9.6 hr. As this is a first order reaction, and the primary
203 limiting step in exchange is lipid desorption, $k_\mathrm{b}$ for \ac{DMPC} is
204 $k_{\mathrm{b}_\mathrm{PC}}=\frac{\log 2}{9.6 \times 60 \times 60} \approx
205 \Sexpr{log(2)/(9.6*60*60)}
206 \mathrm{s}^{-1}$. We assume that $k_\mathrm{b}$ for \ac{SM} is the same as for \ac{PC}.
207 To estimate the $k_\mathrm{b}$ of \ac{PE} and \ac{PS}, we used the data from
208 \citet{Nichols1982:ret_amphiphile_transfer} who measured the rate of
209 exchange of the fluorescent label \ac{C6NBD} attached to different
210 lipid species. Although the values of $k_\mathrm{b}$ are different for the labeled
211 and unlabeled lipids, we assume that the ratios of the kinetics
212 constants for the lipid types are the same. Furthermore we assume that
213 PG behaves similarly to \ac{PS}. Thus, we can determine the $k_\mathrm{b}$ of \ac{PE} and
214 \ac{PS} from the already known $k_\mathrm{b}$ of \ac{PC}. For \ac{C6NBD} labeled \ac{PC},
215 \citet{Nichols1982:ret_amphiphile_transfer} obtained a $k_\mathrm{b}$ of
216 $0.89$~$\mathrm{min}^{-1}$, \ac{PE} of $0.45$~$\mathrm{min}^{-1}$ and PG of
217 $0.55$~$\mathrm{min}^{-1}$. Assuming the ratios of the rate of
218 exchange are the same for unlabeled lipids and labeled lipids, we can
219 determine the $k_\mathrm{b}$ of \ac{PE} and \ac{PS} from the already known $k_\mathrm{b}$ of
220 \ac{PC}~\citep{Wimley1990:dmpc_exchange}. Calculating the ratio, this leads
221 us to $k_{\mathrm{b}_\mathrm{PE}} =
222 \frac{k_{\mathrm{b}_\mathrm{PC}}\times\mathrm{PE}}{\mathrm{PC}} \approx
223 \frac{2\times 10^{-5}\,\mathrm{s}^{-1} \times
224 0.45\,\mathrm{min}^{-1}}{0.89\,\mathrm{min}^{-1}} \approx
225 \Sexpr{log(2)/(9.6*60*60)*0.45/0.89}$~$\mathrm{s}^{-1}$
226 and likewise, $k_{\mathrm{b}_\mathrm{PS}}\approx
227 \Sexpr{log(2)/(9.6*60*60)*0.55/0.89}$~$\mathrm{s}^{-1}$.
228 The $k_\mathrm{b}$ of \ac{SM} was determined using the work of
229 \citet{Bai1997:lipid_movementbodipy}, who measured spontaneous
230 transfer of C$_5$-DMB-SM and C$_5$-DMB-PC from donor and acceptor
231 vesicles, finding $3.4\times10^{-2}$~$\mathrm{s}^{-1}$ and
232 $2.2\times10^{-3}$~$\mathrm{s}^{-1}$ respectively; using the ratio of
233 $k_\mathrm{b}$ of C$_5$-DMB-SM to the $k_\mathrm{b}$ of C$_5$-DMB-PC times the $k_\mathrm{b}$ of
234 \ac{PC} ($\frac{3.4 \times 10^{-2}\mathrm{s}^{-1}}{2.2 \times
235 10^{-3}\mathrm{s}^{-1}}\approx
236 \Sexpr{log(2)/(9.6*60*60)}\mathrm{s}^{-1}$),
237 we obtain $k_{\mathrm{b}_\mathrm{SM}} \approx
238 \Sexpr{log(2)/(9.6*60*60)*3.4e-2/2.2e-3}$.
239 In the case of \ac{CHOL}, \citet{Jones1990:lipid_transfer} measured the
240 $t_{1/2}$ of [$^3$H] transfer from POPC vesicles and found it to be 41
241 minutes, leading to a $k_{\mathrm{b}_\mathrm{CHOL}} = \frac{\log 2}{41\times
243 \Sexpr{log(2)/(41*60)}$~$\mathrm{s}^{-1}$.
245 \subsubsection{Headgroup Surface Area for lipid types}
247 % Smondyrev, A. M., and M. L. Berkowitz. 1999b. United atom force
248 % field for phospholipid membranes: constant pressure molecular
249 % dynamics simulation of dipalmitoylphosphatidicholine/water system.
250 % J. Comput. Chem. 50:531–545.
252 Different lipids have different headgroup surface areas, which contributes to
253 $\left[S_\mathrm{vesicle}\right]$. \citet{Smaby1997:pc_area_with_chol}
254 measured the surface area of POPC with a Langmuir film balance, and
255 found it to be 63~Å$^2$ at $30$~$\frac{\mathrm{mN}}{\mathrm{m}}$.
256 Molecular dynamic simulations found an area of 54 Å$^2$ for
257 DPPS\citep{Cascales1996:mds_dpps_area,Pandit2002:mds_dpps}, which is
258 in agreement with the experimental value of 56~Å$^2$ found using a
259 Langmuir balance by \citet{Demel1987:ps_area}.
260 \citet{Shaikh2002:pe_phase_sm_area} measured the area of \ac{SM} using a
261 Langmuir film balance, and found it to be 61~Å$^2$. Using $^2$H NMR,
262 \citet{Thurmond1991:area_of_pc_pe_2hnmr} found the area of
263 DPPE-d$_{62}$ to be 55.4 Å$^2$. \citet{Robinson1995:mds_chol_area}
264 found an area for \ac{CHOL} of 38~Å$^2$ using molecular dynamic
267 % robinson's chol area is kind of crappy; they did it using MDS, but
268 % vaguely; they don't indicate what it was at the end, and they only
269 % had 2 molecules of \ac{CHOL}.
273 % \ac{PC}: 63 A (Smaby97rd) (Pandit (?)
274 % \ac{PS}: 54 A (Pandit02LIB) (Cascales 1996; J. Chem. Phys 104:2713, Mukhopadhyay2004)
275 % \ac{CHOL}: 38 A (Robinson 1995) (Lundberg 1982)
276 % \ac{SM}: 51 A (Shaikh2002, but 61(?))
277 % PE: 55 A (Thurmond, 1991) (Pandit2002)
279 % Compare to results by Petrache2004 with DOPC of 72.5, DOPS of 65.3.
281 % The cross‑sectional area of cholesterol is ~37 ֵ2 (Lundberg, B. 1982. A
282 % surface film study of the lateral packing of phosphatidylcholine and
283 % cholesterol. Chem. Phys. Lipids. 31:23-32).
285 % From Pandit02LIB: the area per headgroup for PE is ~10‑20\% (Thurmond
286 % et al., 1991) smaller than the area per PC molecule. For DLPE, it is
287 % 50.6 ֵ2 (McIntosh and Simon, 1986; Damodaran et al., 1992). The
288 % estimated area per molecule for DPPE is ~55.4 ֵ2 (Thurmond et al.,
291 % From Pandit02LIB: (from their simulations) average area per headgroup
292 % for DPPS molecules is 53.75 ± 0.10 ֵ2. The values inferred from the
293 % experiments are in the region between 45 and 55 ֵ2 (Cevc et al., 1981.
294 % Titration of the phase transition of phosphatidylserine bilayer
295 % membranes. Effect of pH, surface electrostatics, ion binding, and
296 % head‑group hydration. Biochemistry. 20:4955‑4965; | Demel et al.,
297 % 1987. Monolayer characteristics and thermal behavior of natural and
298 % synthetic phosphatidylserines. Biochemistry. 26:8659‑8665). Note that
299 % the average area per headgroup for DPPC bilayer obtained from
300 % simulations and measurements is ~62 ֵ A2. One would expect that the
301 % area per headgroup in case of DPPS molecules will be larger than the
302 % area per DPPC molecule, because the headgroups of DPPS are charged and
303 % therefore repel each other. Contrary to this expectation, the area per
304 % headgroup in DPPS is ~13\% smaller than that of DPPC. López Cascales
305 % et al. (1996. Molecular dynamics simulation of a charged biological
306 % membrane. J. Chem. Phys. 104:2713‑2720.) proposed that this reduction
307 % in the area per headgroup is due to a strong intermolecular
308 % coordination between DPPS molecules. MD gave area per POPS of 55 A2 at
309 % 300K (Mukhopadhyay04LIB).
311 % However, by 2H NMR and X-ray, Petrache04LIB determined the area of
312 % DOPS at 30 C to be 65.3 A2, considerably smaller than DOPC, determined
313 % by the same group in another paper to be 72.5 ֵ A2
314 % (Tristram‑Nagle00LIB).
316 % From Shaikh02: At 35 mN/m of surface pressure, the areas/molecule for
317 % POPE, SM, and CHOL were found to be 65.2, 61.1, and 38.8 Å2.
319 % In Lonchin99: POPC area 72 A2 ((a) Cornell, B. A.; Middlehurst, J.;
320 % Separovic, F. Biochim. Biophys. Acta1980, 598, 405−410), oleic acid
321 % area 32 A2 ((a) Small, D. M. In Handbook of Lipid Research; Hanahan,
322 % D. J., Ed.; Plenum Press: New York, 1986; Vol. 4, Chapter 9.
324 % (b) Jain, M. K. Introduction to Biological Membranes, 2nd ed.; John
325 % Wiley \& Sons: New York, 1988; Chapter 3). Cistola, D. P.; Atkinson,
326 % D.; Hamilton, J. A.; Small, D. M. Biochemistry1986, 25, 2804−2812.
328 % From Marsh96: The values obtained for the headgroup area per molecule,
329 % fitted as free parameters in Fig. 6, are APE = 58 A2 and APC = 78 A2,
330 % for DOPE and DOPC, respectively. These values are within the range
331 % obtained for the corresponding areas per molecule of
332 % phosphatidylethanolamines and phosphatidylcholines in the fluid
333 % lamellar phase (see, e.g., Marsh, 1990), and correspond to values for
334 % the lipid packing parameter of V/Al = 1.08 and 1.34 for DOPC and DOPE,
335 % respectively. The former value is consistent with DOPC alone forming
336 % hydrated lamellar structures rather than nonlamellar structures. The
337 % latter value lies at the lower end of the packing parameters obtained
338 % for phosphatidylethanolamines.
340 % From Kumar91: The head group area chosen was 71.7 A2 for PC and
341 % lysoPC, 42 A2 for PE, and 19 A2 for Chol (5, 14). The value selected
342 % for PC is in agreement with the 71 A2 value determined for PCs with
343 % saturated acyl chains (14). The hydrocarbon volume and the length
344 % chosen for Chol are 400 A3 and 17.25 A, respectively (5).
348 % From Sampaio05: Besides this work and our own earlier report on the
349 % association of NBD-DMPE with lipid bilayers (Abreu et al., 2004), we
350 % are aware of only one other report in the literature (Nichols, 1985)
351 % in which all the kinetic constants of lipid-derived amphiphile
352 % association with lipid bilayer membranes were experimentally measured.
353 % {indeed; everything is k- !!!; rz}
355 % From McLean84LIB: Although it is difficult to measure cmc values for
356 % the sparingly soluble lipids used in this study, estimates for
357 % lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds,
358 % 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution;
359 % Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6
360 % X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X
361 % 10-8 M for DMPC was estimated from the linear relationship between ln
362 % cmc and the number of carbons in the PC acyl chain by using data for n
363 % = 7, 8, 10, and 16 [summarized in Tanford (1980)].
365 % From Toyota08: Recently, several research groups have reported
366 % self-reproducing systems of giant vesicles that undergo a series of
367 % sequential transformations: autonomous growth, self-division, and
368 % chemical reactions to produce membrane constituents within the giant
371 % Vesicle sizes of 25 nm for SUV and 150 nm for LUV were mentioned by
374 % From Lund-Katz88: Charged and neutral small unilamellar vesicles
375 % composed of either saturated PC, unsaturated PC, or SM had similar
376 % size distributions with diameters of 23 \& 2 nm.
378 % From Sampaio05LIB: The exchange of lipids and lipid derivatives
379 % between lipid bilayer vesicles has been studied for at least the last
380 % 30 years. Most of this work has examined the exchange of amphiphilic
381 % molecules between a donor and an acceptor population. The measured
382 % efflux rates were shown in almost all cases, not surprisingly, to be
383 % first order processes. In all of this work, the insertion rate has
384 % been assumed to be much faster than the efflux rate. Having measured
385 % both the insertion and desorption rate constants for amphiphile
386 % association with membranes, our results show that this assumption is
387 % valid. In several cases reported in the literature, the insertion rate
388 % constant was assumed, although never demonstrated, to be a
389 % diffusion-controlled process.
391 % (for methods? From McLean84LIB: The activation free energies and free
392 % energies of transfer from self-micelles to water increase by 2.2 and
393 % 2.1 kJ mol-' per methylene group, respectively. {see if we can use it
394 % to justify arranging our changed in activating energy around 1
397 % Jones90 give diameter of LUV as 100 nm, and of SUV as 20 nm; that
398 % would give the number of molecules per outer leaflet of a vesicle as
401 % Form Simard08: Parallel studies with SUV and LUV revealed that
402 % although membrane curvature does have a small effect on the absolute
403 % rates of FA transfer between vesicles, the ΔG of membrane desorption
404 % is unchanged, suggesting that the physical chemical properties which
405 % govern FA desorption are dependent on the dissociating molecule rather
406 % than on membrane curvature. However, disagreements on this fundamental
407 % issue continue (57, 61, 63, 64)
409 % (methods regarding the curvature effect: Kleinfeld93 showed that the
410 % transfer parameters of long-chain FFA between the lipid vesicles
411 % depend on vesicle curvature and composition. Transfer of stearic acid
412 % is much slower from LUV as compared to SUV).
414 % From McLean84: In a well-defined experimental system consisting of
415 % unilamellar lipid vesicles, in the absence of protein, the
416 % rate-limiting step for the overall exchange process is desorption
417 % (McLean \& Phillips, 1981). {thus I can take exchange data for the
418 % estimation of k- rz; 8/11/08}.
420 \subsubsection{Complex Formation 1}
422 \citet{Thomas1988:chol_transfer} found that \ac{SM} significantly decreases
423 the rate of cholesterol transfer from \ac{PC}, while \ac{PE} and \ac{PS} have no
424 effect at physiologically significant levels. We attribute this to the
425 formation of complexes between \ac{SM} and \ac{CHOL}, which retards the rate of
426 efflux of \ac{CHOL} from the membrane. Similar complexes exist between \ac{PC}
427 and \ac{CHOL}, where it was shown that two \ac{CHOL} molecules complex with one
428 \ac{PC}~\citep{Huang1999:chol_solubility_model,
429 Huang1999:cholesterol_solubility,McConnell2006,McConnell2003}. It
430 has also been shown that \ac{SM} binds more \ac{CHOL} molecules than
431 \ac{PC}~\citep{Katz1988:pl_packing_chol}. We assume that three \ac{CHOL} complex
432 with one \ac{SM}, leading to $\mathrm{CF}1$ values of $3$, $2$, and $-1$
433 for \ac{SM}, \ac{PC}, and \ac{CHOL}, respectively.
435 \subsubsection{Curvature}
437 We used the formulation for curvature of
438 \citet{Israelachvili1975:amphiphile_self_assembly}, or $S=\frac{v}{l_c
439 a}$, where $S$ is the curvature, which ranges from $0$ to $\infty$,
440 $l_c$ is the critical length of the acyl chain, $v$ is the volume of
441 the lipid, and $a$ is the area of the head group.
442 \citet{Kumar1991:lipid_packing} found a value for $S$ of 1.33 for PE,
443 1.21 for \ac{CHOL}, and 0.8 for diC$_{16}$ \ac{PC}. We assume that \ac{PS} has neutral
444 curvature (value of 1), and that \ac{SM} has the same curvature as \ac{PC}
445 (0.8). It should also be noted that in reality the curvature parameter
446 changes with length, but at longer chain lengths, is effectively
447 constant; in the current model, curvature is held constant for each
451 % 1. Israelachvili, J. N., Mitchell, D. J. & Ninham, B. W. (1976) J.
452 % Chem. Soc. Faraday Trans. 2 72, 1525-1568.
453 % 2. Israelachvili, J. N., Marcelja, S. & Horn, R. G. (1980) Q. Rev.
454 % Biophys. 13, 121-200.
455 % 3. Israelachvili, J. N. (1985) Intermolecular and Surface Forces
456 % (Academic, New York), pp. 229-271.
458 % \DLA{This section needs to be added still.}
460 % PE = 1.33 (Kumar91)
461 % CHOL = 1.21 (Kumar91)
463 % SM=0.8 (assumed by rz same as PC)
464 % PS=1 (no refs so far; should be close to unity; rz)
466 % From Epand94LIB: a considerable amount of heat is evolved per mole of
467 % bilayer stabilizer when such molecules are inserted into membranes
468 % which are under a high curvature strain. If this energy were coupled
469 % to a membrane event, such as the conformational change in a protein,
470 % it could be the driving force responsible for such processes.
472 % From Epand94LIB: considerable energy may be released upon the
473 % incorporation of certain molecules into membranes which have a low
474 % radius of spontaneous curvature. {discuss with DL in terms of
475 % catalytic activity; rz; 8/24/2010}.
477 % Kleinfeld93 showed that the transfer parameters of long-chain FFA
478 % between the lipid vesicles depend on vesicle curvature and
479 % composition. Transfer of stearic acid is much slower from LUV as
482 % Substitution of DPPC SUV for POPC SUV as the donor matrix revealed
483 % some differences in the energetics of transfer that are associated
484 % with the physical state of the donor particle. Whereas the free
485 % energies of transfer of oleic acid from DPPC SUV and POPC SUV are
486 % similar, the free energy of transfer of stearic acid from DPPC SUV is
487 % much higher than that observed from POPC.
489 % From Kumar91: S is given by S = V/al, where V is the hydrocarbon
490 % volume, a is the area of head group, and l is the critical length of
491 % the hydrocarbon chain (1-3). a, V, and I are all estimable or
492 % measurable (4), and the value of S can be calculated. The value of S
493 % determines the aggregate formed by lipids or any amphiphiles upon
494 % hydration. It has been shown that lipids aggregate to form spherical
495 % micelles (S < 1/3), nonspherical (cylindrical) micelles (1/3 < S <
496 % 1/2), bilayers (1/2 < S < 1), and reverse micelles or hexagonal (HII,)
497 % phases (S > 1). However, theoreticians caution that the above
498 % predicted limits set on the values of S are relatively insensitive to
499 % the exact values of V and a but are strongly dependent upon the choice
502 % From Kumar91: Cylindrical and wedge-shaped molecules have been shown
503 % to be essential for spontaneous vesiculation (6) as well as bilayer
504 % stabilization (7, 8).
506 % From Kumar91: It has also been shown experimentally that molecules
507 % having complementary molecular shapes can form bilayer structures.
508 % Inverted cone (wedge)- shaped lysoPC and cone-shaped Chol or
509 % unsaturated PE can interact stoichiometrically to form bilayer
510 % structures (9-12). {quote to justify our formula for complementary
513 \section{Kinetic adjustments}
514 \label{sec:kinetic_adjustments}
516 In the mass action equation used in the formalism, both the forward
517 and backwards kinetic constants ($k_\mathrm{f}$ and $k_\mathrm{b}$) are adjusted (by
518 $k_{\mathrm{f}i\mathrm{adj}}$ and $k_{\mathrm{b}i\mathrm{adj}}$) to reflect the
519 properties of the vesicle and the properties of the monomer entering
520 or exiting the vesicle.
522 \subsection{Forward Kinetic Adjustments}
524 The forward rate constant adjustment, $k_{\mathrm{f}i\mathrm{adj}}$ takes into
525 account unsaturation ($un_\mathrm{f}$), charge ($ch_\mathrm{f}$), curvature ($cu_\mathrm{f}$),
526 length ($l_\mathrm{f}$), and complex formation ($CF1_\mathrm{f}$), each of which is
527 modified depending on the specific component and the vesicle into
528 which the component is entering; the final $k_{\mathrm{f}i\mathrm{adj}}$ is the
529 product of the adjustments made separately for each property.
532 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}
536 \subsubsection{Unsaturation Forward}
538 In order for a lipid to be inserted into a membrane, a void has to be
539 formed for it to fill. Voids can be generated by the combination of
540 unsaturated and saturated lipids forming heterogeneous domains. Void
541 formation is increased when the unsaturation of lipids in the vesicle
542 is widely distributed; in other words, the insertion of lipids into
543 the membrane is greater when the standard deviation of the
544 unsaturation is larger. %
545 %%% \RZ{May need to site (at least for us) works showing
546 %%% mismatch-dependent ``defects''}. %
547 Assuming that an increase in width of the distribution linearly
548 decreases the free energy of activation, the $un_\mathrm{f}$ parameter must
549 follow $a^{\mathrm{stdev}\left(un_\mathrm{vesicle}\right)}$ where $a > 1$,
550 so a convenient starting base for $a$ which leads to a reasonable
551 change in Eyring activation energy is $2$:
554 un_\mathrm{f} = 2^{\mathrm{stdev}\left(un_\mathrm{vesicle}\right)}
555 \label{eq:unsaturation_forward}
558 The mean $\mathrm{stdev}\left(un_\mathrm{vesicle}\right)$ in our
559 simulations is around $1.5$, which leads to a $\Delta \Delta
560 G^\ddagger$ of $\Sexpr{to.latex(format(digits=3,to.kcal(2^1.5)))}
561 \frac{\mathrm{kcal}}{\mathrm{mol}}$, and a total range of possible
562 values depicted in \cref{fig:unf_graph}.
564 % \RZ{Explain here, or even earlier that the formulas were ad hoc
565 % adjusted to correspond to ``reasonable'' changes in the Eyring
566 % activation energies.}
568 % It is not clear that the unsaturation of the inserted monomer will
569 % affect the rate of the insertion positively or negatively, so we do
570 % not include a term for it in this formalism.
574 <<unf_graph,echo=FALSE,results='hide',fig.width=10,fig.height=5,cache=FALSE>>=
575 layout(matrix(1:2,nrow=1,ncol=2))
576 curve(2^x,from=0,to=sd(c(0,4)),
577 xlab=expression(paste(stdev,group("(",italic(un[vesicle]),")"))),
578 ylab=expression(italic(un[f])))
579 vps <- baseViewports()
580 pushViewport(vps$figure)
582 just=c("left","top"),
583 gp=gpar(fontsize=48),
584 y=unit(1,"npc")-unit(0.25,"lines"),
585 x=unit(0,"npc")+unit(0.25,"lines"))
587 curve(to.kcal(2^x),from=0,to=sd(c(0,4)),
588 xlab=expression(paste(stdev,group("(",italic(un[vesicle]),")"))),
589 ylab="Activation Energy (kcal/mol)")
590 vps <- baseViewports()
591 pushViewport(vps$figure)
593 just=c("left","top"),
594 gp=gpar(fontsize=48),
595 y=unit(1,"npc")-unit(0.25,"lines"),
596 x=unit(0,"npc")+unit(0.25,"lines"))
599 \caption{Change in unsaturation forward ($un_\mathrm{f}$) (A) or activation
600 energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the
601 standard deviation of the unsaturation of the vesicle
602 ($\mathrm{stdev}\left(un_\mathrm{vesicle}\right)$).}
603 \label{fig:unf_graph}
606 \subsubsection{Charge Forward}
608 A charged lipid such as \ac{PS} approaching a vesicle with an average
609 charge of the same sign will experience repulsion, whereas those with
610 different signs will experience attraction, the degree of which is
611 dependent upon the charge of the monomer and the average charge of the
612 vesicle. If either the vesicle or the monomer has no charge, there
613 should be no effect of charge upon the rate. This leads us to the
614 following equation, $a^{-\left<ch_v\right> ch_m}$ (which is similar to
615 the numerator of Coulomb's law), where $\left<ch_v\right>$ is the
616 average charge of the vesicle, and $ch_m$ is the charge of the
617 monomer. If either $\left<ch_v\right>$ or $ch_m$ is 0, the adjustment
618 parameter is 1 (no change), whereas it decreases if both are positive
619 or negative, as the product of two real numbers with the same sign is
620 always positive. The base for $a$ ($60$) was chosen ad hoc to
621 correspond to a maximum of about
622 $0.5\frac{\mathrm{kcal}}{\mathrm{mol}}$ change in activation energy
623 for the common case, resulting in the following equation:
627 ch_\mathrm{f} = 60^{-\left<{ch}_v\right> {ch}_m}
628 \label{eq:charge_forward}
631 The most common $\left<{ch}_v\right>$ is around $-0.165$, which leads
632 to a range of $\Delta \Delta G^\ddagger$ from
633 $\Sexpr{format(digits=3,to.kcal(60^(-.165*-1)))}
634 \frac{\mathrm{kcal}}{\mathrm{mol}}$ to
635 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$, and the total range of possible
636 values seen in \cref{fig:chf_graph}.
640 <<chf_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
641 trellis.device(new=F,color=TRUE)
642 trellis.par.set(list(axis.line =list(col="transparent")))
643 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
644 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
645 x <- seq(-1,0,length.out=20)
646 y <- seq(-1,0,length.out=20)
647 grid <- expand.grid(x=x,y=y)
648 grid$z <- as.vector(60^(-outer(x,y)))
649 print(wireframe(z~x*y,grid,cuts=50,
652 scales=list(arrows=F,col=1,distance=0.75),
653 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
654 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
655 zlab=list(label=expression(italic(ch[f])),rot=93)),newpage=FALSE)
658 just=c("left","top"),
659 gp=gpar(fontsize=48),
660 y=unit(1,"npc")-unit(0.25,"lines"),
661 x=unit(0,"npc")+unit(0.25,"lines"))
663 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
664 x <- seq(-1,0,length.out=20)
665 y <- seq(-1,0,length.out=20)
666 grid <- expand.grid(x=x,y=y)
667 grid$z <- as.vector(to.kcal(60^(-outer(x,y))))
668 print(wireframe(z~x*y,grid,cuts=50,
671 scales=list(arrows=F,col=1,distance=0.75),
672 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
673 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
674 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
677 just=c("left","top"),
678 gp=gpar(fontsize=48),
679 y=unit(1,"npc")-unit(0.25,"lines"),
680 x=unit(0,"npc")+unit(0.25,"lines"))
683 \caption{Change in charge forward ($ch_\mathrm{f}$) (A) or activation energy in
684 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the charge of the
685 monomer entering ($ch_\mathrm{monomer}$) and the average charge of
686 the vesicle ($\left<ch_\mathrm{vesicle}\right>$).}
687 \label{fig:chf_graph}
690 \subsubsection{Curvature Forward}
692 Curvature is a measure of the intrinsic propensity of specific lipids
693 to form micelles (positive curvature), inverted micelles (negative
694 curvature), or planar sheets (neutral
695 curvature)~\citep{Israelachvili1975:amphiphile_self_assembly}. %
696 In this formalism, curvature is measured as the ratio of the volume of
697 the lipid to the area of the head times the length of the lipid
698 ($S=\frac{v}{l_c a}$), so negative curvature is bounded by $(0,1)$,
699 neutral curvature is 1, and positive curvature is bounded by
700 $(1,\infty)$. The curvature can be transformed using $\log$, which has
701 the property of making the range of positive and negative curvature
702 equal, and distributed about 0.
704 As in the case of unsaturation, void formation is increased by the
705 presence of lipids with mismatched curvature. Thus, a larger
706 distribution of curvature in the vesicle increases the rate of lipid
707 insertion into the vesicle. However, a component with curvature
708 $cu^{-1}$ will cancel out a component with curvature $cu$, so we have
709 to log transform (turning these into $-\log cu$ and $\log cu$), then
710 take the absolute value ($\log cu$ and $\log cu$), and finally measure
711 the width of the distribution, which in the case of exactly mismatched
712 curvatures is $0$. Thus, by using the log transform to make the range
713 of the lipid curvature equal between positive and negative, and taking
714 the average to cancel out exactly mismatched curvatures, we come to an
715 equation with the shape $a^{\left<\log cu_\mathrm{vesicle}\right>}$.
716 An ad hoc base for $a$ which for the common case is $10$, yielding:
720 % cu_\mathrm{f} = 10^{\mathrm{stdev}\left|\log cu_\mathrm{vesicle}\right|}
721 cu_\mathrm{f} = 10^{\left|\left<\log cu_\mathrm{monomer} \right>\right|\mathrm{stdev} \left|\log cu_\mathrm{vesicle}\right|}
722 \label{eq:curvature_forward}
725 The most common $\left|\left<\log {cu}_v\right>\right|$ is around
726 $0.013$, which with the most common $\mathrm{stdev} \log
727 cu_\mathrm{vesicle}$ of $0.213$ leads to a $\Delta \Delta G^\ddagger$
728 of $\Sexpr{format(digits=3,to.kcal(10^(0.13*0.213)))}
729 \frac{\mathrm{kcal}}{\mathrm{mol}}$. This is a consequence of the
730 relatively matched curvatures in our environment. The full range of
731 $cu_\mathrm{f}$ values possible are shown in \cref{fig:cuf_graph}.
733 % 1.5 to 0.75 3 to 0.33
735 <<cuf_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
736 trellis.device(new=F,color=TRUE)
737 trellis.par.set(list(axis.line =list(col="transparent")))
738 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
739 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
740 grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
741 sd(abs(log(c(1,0.33)))),
742 sd(abs(log(c(0.33,3)))))),length.out=20),
743 y=seq(0,max(c(mean(log(c(1,3)),
744 mean(log(c(1,0.33))),
745 mean(log(c(0.33,3)))))),length.out=20))
746 grid$z <- 10^(grid$x*grid$y)
747 print(wireframe(z~x*y,grid,cuts=50,
750 scales=list(arrows=F,col=1,distance=0.75),
751 xlab=list(label=expression(paste(stdev,group("(",italic(log(cu[vesicle])),")"))),rot=30),
752 ylab=list(label=expression(paste(symbol("\341"),italic(log(cu[vesicle])),symbol("\361"))),rot=-35),
753 zlab=list(label=expression(italic(cu[f])),rot=93)),newpage=FALSE)
756 just=c("left","top"),
757 gp=gpar(fontsize=48),
758 y=unit(1,"npc")-unit(0.25,"lines"),
759 x=unit(0,"npc")+unit(0.25,"lines"))
761 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
762 grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
763 sd(abs(log(c(1,0.33)))),sd(abs(log(c(0.33,3)))))),length.out=20),
764 y=seq(0,max(c(mean(log(c(1,3)),
765 mean(log(c(1,0.33))),
766 mean(log(c(0.33,3)))))),length.out=20))
767 grid$z <- to.kcal(10^(grid$x*grid$y))
768 print(wireframe(z~x*y,grid,cuts=50,
771 scales=list(arrows=F,col=1,distance=0.75),
772 xlab=list(label=expression(paste(stdev,group("(",italic(log(cu[vesicle])),")"))),rot=30),
773 ylab=list(label=expression(paste(symbol("\341"),italic(log(cu[vesicle])),symbol("\361"))),rot=-35),
774 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
777 just=c("left","top"),
778 gp=gpar(fontsize=48),
779 y=unit(1,"npc")-unit(0.25,"lines"),
780 x=unit(0,"npc")+unit(0.25,"lines"))
783 \caption{Change in curvature forward ($cu_\mathrm{f}$) (A) or activation energy
784 in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the standard
785 deviation of the log of the curvature of the vesicle
786 ($\mathrm{stdev}\left(\log cu_\mathrm{vesicle}\right)$) and the mean
787 of the log of the curvature of the vesicle ($\left<\log
788 cu_\mathrm{vesicle}\right>$).}
789 \label{fig:cuf_graph}
792 \subsubsection{Length Forward}
794 As in the case of unsaturation, void formation is easier when vesicles
795 are made up of components of mismatched lengths. Thus, when the width
796 of the distribution of lengths is larger, the forward rate should be
797 greater as well, leading us to an equation of the form
798 $x^{\mathrm{stdev} l_\mathrm{vesicle}}$, where $\mathrm{stdev}
799 l_\mathrm{vesicle}$ is the standard deviation of the length of the
800 components of the vesicle, which has a maximum possible value of
801 $\Sexpr{format(digits=3,sd(c(rep(12,50),rep(24,50))))}$ and a minimum
802 of 0 in this set of simulations. Based on activation energy, a
803 reasonable base for $x$ is 2, leading to:
806 l_\mathrm{f} = 2^{\mathrm{stdev} l_\mathrm{vesicle}}
807 \label{eq:length_forward}
810 The most common $\mathrm{stdev} l_\mathrm{vesicle}$ is around $3.4$, which leads to
811 a $\Delta \Delta G^\ddagger$ of
812 $\Sexpr{format(digits=3,to.kcal(2^(3.4)))}
813 \frac{\mathrm{kcal}}{\mathrm{mol}}$.
815 % While it could be argued that increased length of the monomer could
816 % affect the rate of insertion into the membrane, it's not clear whether
817 % it would increase (by decreasing the number of available hydrogen
818 % bonds, for example) or decrease (by increasing the time taken to fully
819 % insert the acyl chain, for example) the rate of insertion or to what
820 % degree, so we do not take it into account in this formalism.
822 % \DLA{Incorporate McLean84 here}
823 % From McLean84LIB: Although it is difficult to measure cmc values for
824 % the sparingly soluble lipids used in this study, estimates for
825 % lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds,
826 % 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution;
827 % Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6
828 % X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X
829 % 10-8 M for DMPC was estimated from the linear relationship between ln
830 % cmc and the number of carbons in the PC acyl chain by using data for n
831 % = 7, 8, 10, and 16 [summarized in Tanford (1980)].
833 % From Nichols85: The magnitude of the dissociation rate constant
834 % decreases by a factor of approximately 3.2 per carbon increase in acyl
835 % chain length (see Table II here) {take into account for the formula;
838 % From Nichols85: The magnitude of the partition coefficient increases
839 % with acyl chain length [Keq(M-C6-NBD-PC) = (9.8±2.1) X 106 M-1 and Keq
840 % (P-C6-NBD-PC) = (9.4±1.3) x 107 M-1. {take into account for the
841 % formula; rz 8/17/2010}.
843 % From Nichols85: The association rate constant is independent of acyl
844 % chain length. {take into account for the formula; rz 8/17/2010}.
849 <<lf_graph,echo=FALSE,results='hide',fig.width=14,fig.height=5>>=
850 layout(matrix(1:2,nrow=1,ncol=2))
851 curve(2^x,from=0,to=sd(c(12,24)),
852 xlab=expression(paste(stdev,group("(",italic(l[vesicle]),")"))),
853 ylab=expression(italic(l[f])),
856 vps <- baseViewports()
857 pushViewport(vps$figure)
859 just=c("left","top"),
860 gp=gpar(fontsize=48),
861 y=unit(1,"npc")-unit(0.25,"lines"),
862 x=unit(0,"npc")+unit(0.25,"lines"))
864 curve(to.kcal(2^x),from=0,to=sd(c(12,24)),
865 xlab=expression(paste(stdev,group("(",italic(l[vesicle]),")"))),
866 ylab="Activation Energy (kcal/mol)",
868 vps <- baseViewports()
869 pushViewport(vps$figure)
871 just=c("left","top"),
872 gp=gpar(fontsize=48),
873 y=unit(1,"npc")-unit(0.25,"lines"),
874 x=unit(0,"npc")+unit(0.25,"lines"))
877 \caption{Change in length forward ($l_\mathrm{f}$) (A) or activation energy in
878 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the standard
879 deviation of the length of the vesicle
880 ($\mathrm{stdev}\left(l_\mathrm{vesicle}\right)$).}
884 \subsubsection{Complex Formation}
885 There is no contribution of complex formation to the forward reaction
886 rate in the current formalism.
890 \label{eq:complex_formation_forward}
893 \subsection{Backward adjustments ($k_{\mathrm{b}i\mathrm{adj}}$)}
895 Just as the forward rate constant adjustment $k_{\mathrm{f}i\mathrm{adj}}$
896 does, the backwards rate constant adjustment $k_{\mathrm{b}i\mathrm{adj}}$
897 takes into account unsaturation ($un_\mathrm{b}$), charge ($ch_\mathrm{b}$), curvature
898 ($cu_\mathrm{b}$), length ($l_\mathrm{b}$), and complex formation ($CF1_\mathrm{b}$), each of
899 which is modified depending on the specific component and the vesicle
900 from which the component is exiting:
904 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}
908 \subsubsection{Unsaturation Backward}
910 Unsaturation also influences the ability of a lipid molecule to leave
911 a membrane. If a molecule has an unsaturation level which is different
912 from the surrounding membrane, it will be more likely to leave the
913 membrane. The more different the unsaturation level is, the greater
914 the propensity for the lipid molecule to leave. However, a vesicle
915 with some unsaturation (eg. average unsaturation of 2) is more
916 favorable for lipids with more unsaturation (eg. 3) than lipids with
917 equivalently less unsaturation (eg. 1), so the difference in energy
918 between unsaturation is not linear. Therefore, an equation with the
919 shape $x^{\left| y^{-\left< un_\mathrm{vesicle}\right>
920 }-y^{-un_\mathrm{monomer}} \right| }$, where
921 $\left<un_\mathrm{vesicle}\right>$ is the average unsaturation of the
922 vesicle and $un_\mathrm{monomer}$ is the unsaturation of the monomer,
923 allows for increasing the efflux of molecules from membranes where
924 they strongly mismatch, while allowing vesicles with greater
925 unsaturation to tolerate greater unsaturation mismatch between
926 monomers and the vesicle. The bases x and y were chosen ad hoc to
927 produce reasonable Eyring energies of activation over the range of
928 unsaturations expected, leading to:
932 un_\mathrm{b} = 7^{1-\left(20\left(2^{-\left<un_\mathrm{vesicle} \right>} - 2^{-un_\mathrm{monomer}} \right)^2+1\right)^{-1}}
933 \label{eq:unsaturation_backward}
936 The most common $\left<un_\mathrm{vesicle}\right>$ is around $1.7$, which
937 leads to a range of $\Delta \Delta G^\ddagger$ from
938 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-0)^2+1))))}
939 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with 0 unsaturation
941 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-4)^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
942 for monomers with 4 unsaturations. See \cref{fig:unb_graph} for
943 the full range of possible values.
947 <<unb_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
948 trellis.device(new=F,color=TRUE)
949 trellis.par.set(list(axis.line =list(col="transparent")))
950 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
951 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
952 grid <- expand.grid(x=seq(0,4,length.out=20),
953 y=seq(0,4,length.out=20))
954 grid$z <- (7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1)))
955 print(wireframe(z~x*y,grid,cuts=50,
958 scales=list(arrows=F,col=1,distance=0.75),
959 xlab=list(label=expression(paste(symbol("\341"),italic(un[vesicle]),symbol("\361"))),rot=30),
960 ylab=list(label=expression(italic(un[monomer])),rot=-35),
961 zlab=list(label=expression(italic(un[b])),rot=93)),newpage=FALSE)
963 just=c("left","top"),
964 gp=gpar(fontsize=48),
965 y=unit(1,"npc")-unit(0.25,"lines"),
966 x=unit(0,"npc")+unit(0.25,"lines"))
968 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
969 grid <- expand.grid(x=seq(0,4,length.out=20),
970 y=seq(0,4,length.out=20))
971 grid$z <- to.kcal((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="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
981 just=c("left","top"),
982 gp=gpar(fontsize=48),
983 y=unit(1,"npc")-unit(0.25,"lines"),
984 x=unit(0,"npc")+unit(0.25,"lines"))
987 \caption{Change in unsaturation backward ($un_\mathrm{b}$) (A) or activation energy in
988 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the unsaturation of the
989 monomer leaving ($un_\mathrm{monomer}$) and the average unsaturation of
990 the vesicle ($\left<un_\mathrm{vesicle}\right>$).}
991 \label{fig:unb_graph}
994 \subsubsection{Charge Backwards}
995 As in the case of monomers entering a vesicle, opposites attract.
996 Monomers leaving a vesicle leave faster if their charge has the same
997 sign as the average charge vesicle. An equation of the form
998 $ch_\mathrm{b} = a^{\left<ch_v\right> ch_m}$ is then appropriate, and
999 using a base of $a=20$ yields:
1002 ch_\mathrm{b} = 20^{\left<{ch}_v\right> {ch}_m}
1003 \label{eq:charge_backwards}
1006 The most common $\left<ch_v\right>$ is around $-0.164$, which leads to
1007 a range of $\Delta \Delta G^\ddagger$ from
1008 $\Sexpr{format(digits=3,to.kcal(20^(-.164*-1)))}
1009 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with charge $-1$ to
1010 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with charge $0$.
1011 See \cref{fig:chb_graph} for the full range of possible values
1016 <<chb_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1017 trellis.device(new=F,color=TRUE)
1018 trellis.par.set(list(axis.line =list(col="transparent")))
1019 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1020 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1021 x <- seq(-1,0,length.out=20)
1022 y <- seq(-1,0,length.out=20)
1023 grid <- expand.grid(x=x,y=y)
1024 grid$z <- as.vector(20^(outer(x,y)))
1025 print(wireframe(z~x*y,grid,cuts=50,
1028 scales=list(arrows=F,col=1,distance=0.75),
1029 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
1030 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
1031 zlab=list(label=expression(italic(ch[b])),rot=93)),newpage=FALSE)
1034 just=c("left","top"),
1035 gp=gpar(fontsize=48),
1036 y=unit(1,"npc")-unit(0.25,"lines"),
1037 x=unit(0,"npc")+unit(0.25,"lines"))
1039 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1040 x <- seq(-1,0,length.out=20)
1041 y <- seq(-1,0,length.out=20)
1042 grid <- expand.grid(x=x,y=y)
1043 grid$z <- to.kcal(as.vector(20^(outer(x,y))))
1044 print(wireframe(z~x*y,grid,cuts=50,
1047 scales=list(arrows=F,col=1,distance=0.75),
1048 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
1049 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
1050 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1053 just=c("left","top"),
1054 gp=gpar(fontsize=48),
1055 y=unit(1,"npc")-unit(0.25,"lines"),
1056 x=unit(0,"npc")+unit(0.25,"lines"))
1059 \caption{Change in charge backward ($ch_\mathrm{b}$) (A) or activation energy in
1060 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the charge of the
1061 monomer leaving ($ch_\mathrm{monomer}$) and the average charge of
1062 the vesicle ($\left<ch_\mathrm{vesicle}\right>$).}
1063 \label{fig:chb_graph}
1067 \subsubsection{Curvature Backwards}
1069 The less a monomer's intrinsic curvature matches the average curvature
1070 of the vesicle in which it is in, the greater its rate of efflux. If
1071 the difference is 0, $cu_\mathrm{f}$ needs to be one. To map negative and
1072 positive curvature to the same range, we also need take the logarithm.
1073 Positive and negative curvature lipids cancel each other out,
1074 resulting in an average curvature of 0; the average of the log also
1075 has this property. Increasing mismatches in curvature increase the
1076 rate of efflux, but asymptotically. An equation which satisfies these
1077 criteria has the form $cu_\mathrm{f} = a^{1-\left(b\left( \left< \log
1078 cu_\mathrm{vesicle} \right> -\log
1079 cu_\mathrm{monomer}\right)^2+1\right)^{-1}}$. An alternative
1080 form would use the absolute value of the difference, however, this
1081 yields a cusp and sharp increase about the curvature equilibrium. We
1082 have chosen bases of $a=7$ and $b=20$ ad hoc, based on activation
1083 energy considerations.
1086 cu_\mathrm{b} = 7^{1-\left(20\left(\left< \log cu_\mathrm{vesicle} \right> -\log cu_\mathrm{monomer} \right)^2+1\right)^{-1}}
1087 \label{eq:curvature_backwards}
1090 The most common $\left<\log cu_\mathrm{vesicle}\right>$ is around
1091 $-0.013$, which leads to a range of $\Delta \Delta G^\ddagger$ from
1092 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(0.8))^2+1))))}
1093 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature 0.8 to
1094 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(1.3))^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
1095 for monomers with curvature 1.3 to
1096 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature near
1097 1. The full range of values possible for $cu_\mathrm{b}$ are shown in
1098 \cref{fig:cub_graph}
1100 % \RZ{What about the opposite curvatures that actually do fit to each
1104 % From Kumar91: It has also been shown experimentally that molecules
1105 % having complementary molecular shapes can form bilayer structures.
1106 % Inverted cone (wedge)- shaped lysoPC and cone-shaped Chol or
1107 % unsaturated PE can interact stoichiometrically to form bilayer
1108 % structures (9-12). {quote to justify our formula for complementary
1113 <<cub_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1114 trellis.device(new=F,color=TRUE)
1115 trellis.par.set(list(axis.line =list(col="transparent")))
1116 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1117 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1118 grid <- expand.grid(x=log(seq(0.8,1.33,length.out=20)),
1119 y=seq(0.8,1.33,length.out=20))
1120 grid$z <- 7^(1-1/(20*(grid$x-log(grid$y))^2+1))
1121 print(wireframe(z~x*y,grid,cuts=50,
1124 scales=list(arrows=F,col=1,distance=0.75),
1125 xlab=list(label=expression(paste(symbol("\341"),log(italic(cu[vesicle])),symbol("\361"))),rot=30),
1126 ylab=list(label=expression(italic(cu[monomer])),rot=-35),
1127 zlab=list(label=expression(italic(cu[b])),rot=93)),newpage=FALSE)
1130 just=c("left","top"),
1131 gp=gpar(fontsize=48),
1132 y=unit(1,"npc")-unit(0.25,"lines"),
1133 x=unit(0,"npc")+unit(0.25,"lines"))
1135 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1136 grid <- expand.grid(x=log(seq(0.8,1.33,length.out=20)),
1137 y=seq(0.8,1.33,length.out=20))
1138 grid$z <- to.kcal(7^(1-1/(20*(grid$x-log(grid$y))^2+1)))
1139 print(wireframe(z~x*y,grid,cuts=50,
1142 scales=list(arrows=F,col=1,distance=0.75),
1143 xlab=list(label=expression(paste(symbol("\341"),log(italic(cu[vesicle])),symbol("\361"))),rot=30),
1144 ylab=list(label=expression(italic(cu[monomer])),rot=-35),
1145 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1148 just=c("left","top"),
1149 gp=gpar(fontsize=48),
1150 y=unit(1,"npc")-unit(0.25,"lines"),
1151 x=unit(0,"npc")+unit(0.25,"lines"))
1154 \caption{Change in curvature backward ($cu_\mathrm{b}$) (A) or activation
1155 energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the
1156 curvature of the monomer leaving ($ch_\mathrm{monomer}$) and the
1157 average of the log of the curvature of the vesicle
1158 ($\left<\log cu_\mathrm{vesicle}\right>$).}
1159 \label{fig:cub_graph}
1163 \subsubsection{Length Backwards}
1165 In a model membrane, the dissociation constant increases by a factor
1166 of approximately 3.2 per carbon decrease in acyl chain length~%
1167 \citep{Nichols1985:phospholipid_monomer_vesicle_thermodynamics}.
1168 Unfortunately, the experimental data known to us only measures the
1169 effect of chain lengths less than or equal to the bulk lipid, not the
1170 effect of lengths exceeding it, and is only known for one bulk lipid
1171 component (DOPC). We assume then, that the increase is in relationship
1172 to the average vesicle, and that lipids with larger acyl chain length
1173 will also show an increase in their dissociation constant.
1176 l_\mathrm{b} = 3.2^{\left|\left<l_\mathrm{vesicle}\right>-l_\mathrm{monomer}\right|}
1177 \label{eq:length_backward}
1180 The most common $\left<l_\mathrm{vesicle}\right>$ is around $17.75$,
1181 which leads to a range of $\Delta \Delta G^\ddagger$ from
1182 $\Sexpr{format(digits=3,to.kcal(3.2^abs(12-17.75)))}
1183 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with length 12 to
1184 $\Sexpr{format(digits=3,to.kcal(3.2^abs(24-17.75)))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
1185 for monomers with length 24 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$
1186 for monomers with length near 18. The full range of possible values of
1187 $l_\mathrm{b}$ are shown in \cref{fig:lb_graph}
1189 % (for methods? From McLean84LIB: The activation free energies and free
1190 % energies of transfer from self-micelles to water increase by 2.2 and
1191 % 2.1 kJ mol-' per methylene group, respectively. \RZ{see if we can use it
1192 % to justify arranging our changed in activating energy around 1
1197 <<lb_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1198 trellis.device(new=F,color=TRUE)
1199 trellis.par.set(list(axis.line =list(col="transparent")))
1200 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1201 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1202 grid <- expand.grid(x=seq(12,24,length.out=20),
1203 y=seq(12,24,length.out=20))
1204 grid$z <- 3.2^(abs(grid$x-grid$y))/1e6
1205 print(wireframe(z~x*y,grid,cuts=50,
1208 scales=list(arrows=F,col=1,distance=0.75),
1209 xlab=list(label=expression(paste(symbol("\341"),italic(l[vesicle]),symbol("\361"))),rot=30),
1210 ylab=list(label=expression(italic(l[monomer])),rot=-35),
1211 zlab=list(label=expression(italic(l[b])%*%10^-6),rot=93)),newpage=FALSE)
1214 just=c("left","top"),
1215 gp=gpar(fontsize=48),
1216 y=unit(1,"npc")-unit(0.25,"lines"),
1217 x=unit(0,"npc")+unit(0.25,"lines"))
1219 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1220 grid <- expand.grid(x=seq(12,24,length.out=20),
1221 y=seq(12,24,length.out=20))
1222 grid$z <- to.kcal(3.2^(abs(grid$x-grid$y)))
1223 print(wireframe(z~x*y,grid,cuts=50,
1226 scales=list(arrows=F,col=1,distance=0.75),
1227 xlab=list(label=expression(paste(symbol("\341"),italic(l[vesicle]),symbol("\361"))),rot=30),
1228 ylab=list(label=expression(italic(l[monomer])),rot=-35),
1229 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1232 just=c("left","top"),
1233 gp=gpar(fontsize=48),
1234 y=unit(1,"npc")-unit(0.25,"lines"),
1235 x=unit(0,"npc")+unit(0.25,"lines"))
1238 \caption{Change in length backwards ($l_\mathrm{b}\times 10^{-6}$) (A) or activation energy
1239 in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the length of the
1240 monomer leaving ($l_\mathrm{monomer}$) and the average length of
1241 lipids in the vesicle ($\left<l_\mathrm{vesicle}\right>$).}
1242 \label{fig:lb_graph}
1246 \subsubsection{Complex Formation Backward}
1248 Complex formation ($CF1$) describes the interaction between \ac{CHOL} and
1249 \ac{PC} or \ac{SM}, where \ac{PC} or \ac{SM} protects the hydroxyl group of \ac{CHOL} from
1250 interactions with water %
1251 \citep{Huang1999:chol_solubility_model,
1252 Huang1999:cholesterol_solubility, McConnell2006, McConnell2003}%
1253 . \ac{PC} ($CF1=2$) can interact with two \ac{CHOL}, and \ac{SM} ($CF1=3$) with three
1254 \ac{CHOL} ($CF1=-1$). If the average of $CF1$ is positive (excess of \ac{SM} and
1255 \ac{PC} with regards to complex formation), components with negative $CF1$
1256 (\ac{CHOL}) will be retained. If average $CF1$ is negative, components with
1257 positive $CF1$ are retained. An equation which has this property is
1258 $CF1_\mathrm{b}=a^{\left<CF1_\mathrm{vesicle}\right>
1259 CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{vesicle}\right>
1260 CF1_\mathrm{monomer}\right|}$, where difference of the exponent is
1261 zero if the average $CF1$ and the $CF1$ of the component have the same
1262 sign, or double the product if the signs are different. Based on
1263 activation energy considerations, we took an ad hoc base for $a$ as
1268 CF1_\mathrm{b}=1.5^{\left<CF1_\mathrm{vesicle}\right> CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{vesicle}\right> CF1_\mathrm{monomer}\right|}
1269 \label{eq:complex_formation_backward}
1272 The most common $\left<CF1_\mathrm{vesicle}\right>$ is around $0.925$,
1273 which leads to a range of $\Delta \Delta G^\ddagger$ from
1274 $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*-1-abs(0.925*-1))))}
1275 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
1277 $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*2-abs(0.925*2))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
1278 for monomers with complex formation $2$ to
1279 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
1280 formation $0$. The full range of possible values for $CF1_\mathrm{b}$ are
1281 depicted in \cref{fig:cf1b_graph}.
1286 <<cf1b_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1287 trellis.device(new=F,color=TRUE)
1288 trellis.par.set(list(axis.line =list(col="transparent")))
1289 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1290 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1291 grid <- expand.grid(x=seq(-1,3,length.out=20),
1292 y=seq(-1,3,length.out=20))
1293 grid$z <- 1.5^(grid$x*grid$y-abs(grid$x*grid$y))
1294 print(wireframe(z~x*y,grid,cuts=50,
1297 scales=list(arrows=F,col=1,distance=0.75),
1298 xlab=list(label=expression(paste(symbol("\341"),italic(CF1[vesicle]),symbol("\361"))),rot=30),
1299 ylab=list(label=expression(italic(CF1[monomer])),rot=-35),
1300 zlab=list(label=expression(italic(CF1[b])),rot=93)),newpage=FALSE)
1303 just=c("left","top"),
1304 gp=gpar(fontsize=48),
1305 y=unit(1,"npc")-unit(0.25,"lines"),
1306 x=unit(0,"npc")+unit(0.25,"lines"))
1308 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1309 grid <- expand.grid(x=seq(-1,3,length.out=20),
1310 y=seq(-1,3,length.out=20))
1311 grid$z <- to.kcal(1.5^(grid$x*grid$y-abs(grid$x*grid$y)))
1312 print(wireframe(z~x*y,grid,cuts=50,
1315 scales=list(arrows=F,col=1,distance=0.75),
1316 xlab=list(label=expression(paste(symbol("\341"),italic(CF1[vesicle]),symbol("\361"))),rot=30),
1317 ylab=list(label=expression(italic(CF1[monomer])),rot=-35),
1318 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1321 just=c("left","top"),
1322 gp=gpar(fontsize=48),
1323 y=unit(1,"npc")-unit(0.25,"lines"),
1324 x=unit(0,"npc")+unit(0.25,"lines"))
1327 \caption{Change in complex formation 1 backwards ($CF1_\mathrm{b}$) (A) or
1328 activation energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus
1329 the complex formation of the monomer leaving
1330 ($CF1_\mathrm{monomer}$) and the average complex formation of the
1331 vesicle ($\left<CF1_\mathrm{vesicle}\right>$).}
1332 \label{fig:cf1b_graph}
1335 \section{Simulation Methodology}
1337 \subsection{Overall Architecture}
1339 The simulations are currently run using a program written in perl with
1340 various modules to handle the subsidiary parts. It produces output for
1341 each generation, including the step immediately preceding and
1342 immediately following a vesicle split (and optionally, each step) that
1343 is written to a state file which contains the state of the vesicle,
1344 environment, kinetic parameters, program invocation options, time, and
1345 various other parameters necessary to recreate the state vector at a
1346 given time. This output file is then read by a separate program in
1347 perl to produce different output which is generated out-of-band;
1348 output which includes graphs and statistical analysis is performed
1349 using R~\citep{R:maincite} (and various grid graphics modules) called
1350 from the perl program.
1352 The separation of simulation and output generation allows refining
1353 output, and simulations performed on different versions of the
1354 underlying code to be compared using the same analysis methods and
1355 code. It also allows later simulations to be restarted from a specific
1356 generation, utilizing the same environment.
1358 Random variables of different distributions are calculated using
1359 Math::Random (0.71), which is seeded for each run using entropy from
1360 the linux kernel's urandom device.
1362 Code is available upon request.
1364 \subsection{Environment Creation}
1366 \subsubsection{Components}
1367 The environment contains concentrations of components. In the current
1368 set of simulations, there are \Sexpr{1+4*5*7} different possible
1369 components, consisting of \ac{PC}, PE, \ac{PS}, \ac{SM}, and \ac{CHOL}, with all lipids
1370 except for \ac{CHOL} having 5 possible unsaturations ranging from 0 to 4,
1371 and 7 lengths from $12,14,...,22$ ($7\cdot
1372 5\cdot4+1=\Sexpr{1+4*5*7}$). In cases where the environment has less
1373 than the maximum number of components, the components are selected in
1374 order without replacement from a randomly shuffled deck of components
1375 (with the exception of \ac{CHOL}) represented once until the desired number
1376 of unique components is obtained. \ac{CHOL} is overrepresented
1377 $\Sexpr{5*7}$ times to be at the level of other lipid types, ensuring
1378 that the probability of \ac{CHOL} being absent in the environment is the
1379 same as the probability of any one of the other lipid types (\ac{PS}, PE,
1380 etc.) being absent. This reduces the number of simulations with a
1381 small number of components which are devoid of \ac{CHOL}.
1383 \subsubsection{Concentrations}
1384 Once the components of the environment have been selected, their
1385 concentrations are determined. In experiments where the environmental
1386 concentration is the same across all lipid components, the
1387 concentration is set to $10^{-10}\mathrm{M}$. In other cases, the
1388 environmental concentration is set to a random number from a gamma
1389 distribution with shape parameter 2 and an average of
1390 $10^{-10}\mathrm{M}$. The base concentration ($10^{-10}\mathrm{M}$)
1391 can also be altered at the initialization of the experiment to
1392 specific values for specific lipid types or components.
1394 The environment is a volume which is the maximum number of vesicles
1395 from a single simulation (4096) times the maximum size of the vesicle
1396 (usually 10000) divided by Avagadro's number divided by the total
1397 environmental lipid concentration, or usually
1398 \Sexpr{4096*10000/6.022E23/141E-10}~L.
1400 \subsection{Initial Vesicle Creation}
1402 Initial vesicles are seeded by selecting lipid molecules from the
1403 environment until the vesicle reaches a specific starting size. The
1404 vesicle starting size has gamma distribution with shape parameter 2
1405 and a mean of the per-simulation specified starting size, with a
1406 minimum of 5 lipid molecules. Lipid molecules are then selected to be
1407 added to the lipid membrane according to four different methods. In
1408 the constant method, molecules are added in direct proportion to their
1409 concentration in the environment. The uniform method adds molecules in
1410 proportion to their concentration in the environment scaled by a
1411 uniform random value, whereas the random method adds molecules in
1412 proportion to a uniform random value. The final method is a binomial
1413 method, which adjust the probability of adding a molecule of a
1414 specific component by the concentration of that component, and then
1415 adds the molecules one by one to the membrane. This final method is
1416 also used in the first three methods to add any missing molecules to
1417 the starting vesicle which were unallocated due to the requirement to
1418 add integer numbers of molecules. (For example, if a vesicle was to
1419 have 10 molecules split evenly between three lipid types, the 10th
1420 molecule would be assigned by randomly choosing between the three
1421 lipid types, weighted by their concentration in the environment.)
1423 \subsection{Simulation Step}
1425 Once the environment has been created and the initial vesicle has been
1426 formed, molecules join and leave the vesicle based on the kinetic
1427 parameters and state equation discussed above until the vesicle splits
1428 forming two progeny vesicles. The program then follows both vesicles
1429 and their progeny until the simulation reaches either the maximum
1430 number of vesicles (usually 4096), or the maximum simulation time
1433 \subsubsection{Calculation of Vesicle Properties}
1435 \label{sec:ves_prop}
1436 $S_\mathrm{vesicle}$ is the surface area of the vesicle, and is the sum of
1437 the surface area of all of the individual lipid components; each lipid
1438 type has a different surface area; we assume that the lipid packing is
1439 optimal, and there is no empty space.
1441 For the other vesicle properties (length, unsaturation, charge, and
1442 curvature), we calculate the mean, the standard deviation, the mean of
1443 the log, the mean of the absolute value of the log, the standard
1444 deviation of the log, and the standard deviation of the absolute value
1445 of the log. All cases, when we take the log, we take the log of the
1446 absolute value, and map $\log(0)$ to $0$. For the purpose of plotting
1447 vesicle properties, we only plot the mean of the property.
1449 \subsubsection{Joining and Leaving of Lipid Molecules}
1451 Determining the number of molecules to add to the lipid membrane
1452 ($n_i$) requires knowing $k_{\mathrm{f}i_\mathrm{adj}}$, the surface area of the
1453 vesicle $S_\mathrm{vesicle}$ (see \cref{sec:ves_prop}), the time interval
1454 $dt$ during which lipids are added, the base $k_{\mathrm{f}i}$, and the
1455 concentration of the monomer in the environment
1456 $\left[C_{i\mathrm{vesicle}}\right]$ (see \cref{eq:state}).
1457 $k_{\mathrm{f}i\mathrm{adj}}$ is calculated (see \cref{eq:kf_adj}) based on the
1458 vesicle properties and their hypothesized effect on the rate (in as
1459 many cases as possible, experimentally based)
1460 (see \cref{sec:kinetic_adjustments} for details). $dt$ can be varied
1461 (see \cref{sec:step_duration}), but for a given step is constant. This
1462 leads to the following:
1464 $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$
1466 In the cases where $n_i > 1$, the integer number of molecules is
1467 added. Fractional $n_i$ or the fractional remainder after the addition
1468 of the integer molecules are the probability of adding a specific
1469 molecule, and are compared to a uniformly distributed random value
1470 between 0 and 1. If the random value is less than or equal to the
1471 fractional part of $n_i$, an additional molecule is added.
1473 Molecules leaving the vesicle are handled in a similar manner, with
1475 $n_i = k_{\mathrm{b}i}k_{\mathrm{b}i_\mathrm{adj}}C_{i_\mathrm{vesicle}}\mathrm{N_A}dt$.
1477 Before addition or removal of molecules, the properties of the vesicle
1478 are calculated; this is done so that molecules can be considered to be
1479 added and removed at the same instant, even though additions are
1480 handled first programmatically. This also avoids cases where a removal
1481 would have resulted in a negative number of molecules for a particular
1484 \subsubsection{Step duration}
1485 \label{sec:step_duration}
1487 $dt$ is the time taken for each step of joining, leaving, and checking
1488 split. In the current implementation, it starts with a value of
1489 $10^{-6}$~$\mathrm{s}$ but this is modified in between each step if the
1490 number of molecules joining or leaving is too large or small. If more
1491 than half of the starting vesicle size molecules join or leave in a
1492 single step, $dt$ is reduced by half. If less than the starting
1493 vesicle size molecules join or leave in 100 steps, $dt$ is doubled.
1494 This is necessary to curtail run times and to automatically adjust to
1495 different experimental runs.
1497 \subsubsection{Vesicle splitting}
1499 If a vesicle has grown to a size which is more than double the
1500 starting vesicle size, the vesicle splits. Molecules are assigned to
1501 the progeny vesicles at random, with each progeny vesicle having an
1502 equal chance of getting a single molecule. The number of molecules to
1503 assign to the first vesicle has binomial distribution with a
1504 probability of an event in each trial of 0.5 and a number of trials
1505 equal to the number of molecules.
1509 The environment, initial vesicle, and the state of the vesicle
1510 immediately before and immediately after splitting are stored
1511 to produce later output.
1513 % \section{Analyzing output}
1515 % The analysis of output is handled by a separate perl program which
1516 % shares many common modules with the simulation program. Current output
1517 % includes simulation progress, summary tables, summary statistics, and
1520 % \subsection{PCA plots}
1522 % Two major groups of axes that contribute to vesicle variation between
1523 % subsequent generations are the components and properties of vesicles.
1524 % Each component in a vesicle is an axis on its own; it can be measured
1525 % either as an absolute number of molecules in each component, or the
1526 % fraction of molecules of that component over the total number of
1527 % molecules; the second approach is often more convenient, as it allows
1528 % vesicles with different number of molecules to be directly compared
1529 % (though it hides any effect of vesicle size).
1531 % In order to visualize the variation between and transition of
1532 % subsequent generations of vesicles from their initial state in the
1533 % simulation, to their final state at the termination of the simulation,
1534 % we plot the projection of the generations onto the two principle PCA
1535 % axes (and alternatively, any pairing of the axes).
1537 % \subsection{Carpet plots}
1539 % Carpet plots show the distance/similarity between the composition or
1540 % properties of all generations in a single run. The difference between
1541 % each group of vesicle can be calculated using Euclidean distance
1542 % (properties and compositions) or H similarity (composition only). We
1543 % must use Euclidean distance for properties because the H distance
1544 % metric is invalid when comparing negative vector elements to positive
1547 % In addition to showing distance or similarity, carpet plots also
1548 % depict propersomes and composomes as square boxes on the diagonals and
1549 % propertypes and compotypes as rectangles off the diagonals, each
1550 % propertype or compotype with a distinct color.
1552 % \subsection{Propersomes, propertypes, composomes and compotypes}
1554 % A generation is considered to be a propersome if it is less than $0.6$
1555 % units (by Euclidean distance of normalized properties) away from the
1556 % generation immediately following and preceding. Likewise, a generation
1557 % is in a composome if its H similarity is more than $0.9$ (by the
1558 % normalized H metric) from the preceding and following generations.
1559 % Propersomes and composomes are then assigned to propertypes and
1560 % compotypes using Paritioning Around Medioids
1561 % (PAM). All values of $k$ between 2 and 15
1562 % (or the number of propersomes and composomes, if that is less than 15)
1563 % are tried, and the clustering with the smallest
1564 % silhouette~\citep{Rousseeuw1987:silhouettes} is chosen as the ideal
1565 % clustering~\citep{Shenhav2005:pgard}.
1568 %\bibliographystyle{unsrtnat}
1569 %\bibliography{references.bib}