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 composition 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 range of $\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 are 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, monomers leaving a
996 vesicle leave faster if their charge has the same sign as the average
997 charge vesicle. An equation of the form $ch_\mathrm{b} = a^{\left<ch_v\right>
998 ch_m}$ is then appropriate, and using a base of $a=20$ yields:
1001 ch_\mathrm{b} = 20^{\left<{ch}_v\right> {ch}_m}
1002 \label{eq:charge_backwards}
1005 The most common $\left<ch_v\right>$ is around $-0.164$, which leads to
1006 a range of $\Delta \Delta G^\ddagger$ from
1007 $\Sexpr{format(digits=3,to.kcal(20^(-.164*-1)))}
1008 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with charge $-1$ to
1009 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with charge $0$.
1010 See \cref{fig:chb_graph} for the full range of possible values
1015 <<chb_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1016 trellis.device(new=F,color=TRUE)
1017 trellis.par.set(list(axis.line =list(col="transparent")))
1018 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1019 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1020 x <- seq(-1,0,length.out=20)
1021 y <- seq(-1,0,length.out=20)
1022 grid <- expand.grid(x=x,y=y)
1023 grid$z <- as.vector(20^(outer(x,y)))
1024 print(wireframe(z~x*y,grid,cuts=50,
1027 scales=list(arrows=F,col=1,distance=0.75),
1028 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
1029 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
1030 zlab=list(label=expression(italic(ch[b])),rot=93)),newpage=FALSE)
1033 just=c("left","top"),
1034 gp=gpar(fontsize=48),
1035 y=unit(1,"npc")-unit(0.25,"lines"),
1036 x=unit(0,"npc")+unit(0.25,"lines"))
1038 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1039 x <- seq(-1,0,length.out=20)
1040 y <- seq(-1,0,length.out=20)
1041 grid <- expand.grid(x=x,y=y)
1042 grid$z <- to.kcal(as.vector(20^(outer(x,y))))
1043 print(wireframe(z~x*y,grid,cuts=50,
1046 scales=list(arrows=F,col=1,distance=0.75),
1047 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
1048 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
1049 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1052 just=c("left","top"),
1053 gp=gpar(fontsize=48),
1054 y=unit(1,"npc")-unit(0.25,"lines"),
1055 x=unit(0,"npc")+unit(0.25,"lines"))
1058 \caption{Change in charge backward ($ch_\mathrm{b}$) (A) or activation energy in
1059 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the charge of the
1060 monomer leaving ($ch_\mathrm{monomer}$) and the average charge of
1061 the vesicle ($\left<ch_\mathrm{vesicle}\right>$).}
1062 \label{fig:chb_graph}
1066 \subsubsection{Curvature Backwards}
1068 The less a monomer's intrinsic curvature matches the average curvature
1069 of the vesicle in which it is in, the greater its rate of efflux. If
1070 the difference is 0, $cu_\mathrm{f}$ needs to be one. To map negative and
1071 positive curvature to the same range, we also need take the logarithm.
1072 Positive and negative curvature lipids cancel each other out,
1073 resulting in an average curvature of 0; the average of the log also
1074 has this property. Increasing mismatches in curvature increase the
1075 rate of efflux, but asymptotically. An equation which satisfies these
1076 criteria has the form $cu_\mathrm{f} = a^{1-\left(b\left( \left< \log
1077 cu_\mathrm{vesicle} \right> -\log
1078 cu_\mathrm{monomer}\right)^2+1\right)^{-1}}$. An alternative
1079 form would use the absolute value of the difference, however, this
1080 yields a cusp and sharp increase about the curvature equilibrium. We
1081 have chosen bases of $a=7$ and $b=20$ ad hoc, based on activation
1082 energy considerations.
1085 cu_\mathrm{b} = 7^{1-\left(20\left(\left< \log cu_\mathrm{vesicle} \right> -\log cu_\mathrm{monomer} \right)^2+1\right)^{-1}}
1086 \label{eq:curvature_backwards}
1089 The most common $\left<\log cu_\mathrm{vesicle}\right>$ is around
1090 $-0.013$, which leads to a range of $\Delta \Delta G^\ddagger$ from
1091 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(0.8))^2+1))))}
1092 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature 0.8 to
1093 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(1.3))^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
1094 for monomers with curvature 1.3 to
1095 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature near
1096 1. The full range of values possible for $cu_\mathrm{b}$ are shown in
1097 \cref{fig:cub_graph}
1099 % \RZ{What about the opposite curvatures that actually do fit to each
1103 % From Kumar91: It has also been shown experimentally that molecules
1104 % having complementary molecular shapes can form bilayer structures.
1105 % Inverted cone (wedge)- shaped lysoPC and cone-shaped Chol or
1106 % unsaturated PE can interact stoichiometrically to form bilayer
1107 % structures (9-12). {quote to justify our formula for complementary
1112 <<cub_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1113 trellis.device(new=F,color=TRUE)
1114 trellis.par.set(list(axis.line =list(col="transparent")))
1115 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1116 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1117 grid <- expand.grid(x=log(seq(0.8,1.33,length.out=20)),
1118 y=seq(0.8,1.33,length.out=20))
1119 grid$z <- 7^(1-1/(20*(grid$x-log(grid$y))^2+1))
1120 print(wireframe(z~x*y,grid,cuts=50,
1123 scales=list(arrows=F,col=1,distance=0.75),
1124 xlab=list(label=expression(paste(symbol("\341"),log(italic(cu[vesicle])),symbol("\361"))),rot=30),
1125 ylab=list(label=expression(italic(cu[monomer])),rot=-35),
1126 zlab=list(label=expression(italic(cu[b])),rot=93)),newpage=FALSE)
1129 just=c("left","top"),
1130 gp=gpar(fontsize=48),
1131 y=unit(1,"npc")-unit(0.25,"lines"),
1132 x=unit(0,"npc")+unit(0.25,"lines"))
1134 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1135 grid <- expand.grid(x=log(seq(0.8,1.33,length.out=20)),
1136 y=seq(0.8,1.33,length.out=20))
1137 grid$z <- to.kcal(7^(1-1/(20*(grid$x-log(grid$y))^2+1)))
1138 print(wireframe(z~x*y,grid,cuts=50,
1141 scales=list(arrows=F,col=1,distance=0.75),
1142 xlab=list(label=expression(paste(symbol("\341"),log(italic(cu[vesicle])),symbol("\361"))),rot=30),
1143 ylab=list(label=expression(italic(cu[monomer])),rot=-35),
1144 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1147 just=c("left","top"),
1148 gp=gpar(fontsize=48),
1149 y=unit(1,"npc")-unit(0.25,"lines"),
1150 x=unit(0,"npc")+unit(0.25,"lines"))
1153 \caption{Change in curvature backward ($cu_\mathrm{b}$) (A) or activation
1154 energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the
1155 curvature of the monomer leaving ($ch_\mathrm{monomer}$) and the
1156 average of the log of the curvature of the vesicle
1157 ($\left<\log cu_\mathrm{vesicle}\right>$).}
1158 \label{fig:cub_graph}
1162 \subsubsection{Length Backwards}
1164 In a model membrane, the dissociation constant increases by a factor
1165 of approximately 3.2 per carbon decrease in acyl chain length~%
1166 \citep{Nichols1985:phospholipid_monomer_vesicle_thermodynamics}.
1167 Unfortunately, the experimental data known to us only measures the
1168 effect of chain lengths less than or equal to the bulk lipid, not the
1169 effect of lengths exceeding it, and is only known for one bulk lipid
1170 component (DOPC). We assume then, that the increase is in relationship
1171 to the average vesicle, and that lipids with larger acyl chain length
1172 will also show an increase in their dissociation constant.
1175 l_\mathrm{b} = 3.2^{\left|\left<l_\mathrm{vesicle}\right>-l_\mathrm{monomer}\right|}
1176 \label{eq:length_backward}
1179 The most common $\left<l_\mathrm{vesicle}\right>$ is around $17.75$,
1180 which leads to a range of $\Delta \Delta G^\ddagger$ from
1181 $\Sexpr{format(digits=3,to.kcal(3.2^abs(12-17.75)))}
1182 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with length 12 to
1183 $\Sexpr{format(digits=3,to.kcal(3.2^abs(24-17.75)))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
1184 for monomers with length 24 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$
1185 for monomers with length near 18. The full range of possible values of
1186 $l_\mathrm{b}$ are shown in \cref{fig:lb_graph}
1188 % (for methods? From McLean84LIB: The activation free energies and free
1189 % energies of transfer from self-micelles to water increase by 2.2 and
1190 % 2.1 kJ mol-' per methylene group, respectively. \RZ{see if we can use it
1191 % to justify arranging our changed in activating energy around 1
1196 <<lb_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1197 trellis.device(new=F,color=TRUE)
1198 trellis.par.set(list(axis.line =list(col="transparent")))
1199 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1200 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1201 grid <- expand.grid(x=seq(12,24,length.out=20),
1202 y=seq(12,24,length.out=20))
1203 grid$z <- 3.2^(abs(grid$x-grid$y))/1e6
1204 print(wireframe(z~x*y,grid,cuts=50,
1207 scales=list(arrows=F,col=1,distance=0.75),
1208 xlab=list(label=expression(paste(symbol("\341"),italic(l[vesicle]),symbol("\361"))),rot=30),
1209 ylab=list(label=expression(italic(l[monomer])),rot=-35),
1210 zlab=list(label=expression(italic(l[b])%*%10^-6),rot=93)),newpage=FALSE)
1213 just=c("left","top"),
1214 gp=gpar(fontsize=48),
1215 y=unit(1,"npc")-unit(0.25,"lines"),
1216 x=unit(0,"npc")+unit(0.25,"lines"))
1218 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1219 grid <- expand.grid(x=seq(12,24,length.out=20),
1220 y=seq(12,24,length.out=20))
1221 grid$z <- to.kcal(3.2^(abs(grid$x-grid$y)))
1222 print(wireframe(z~x*y,grid,cuts=50,
1225 scales=list(arrows=F,col=1,distance=0.75),
1226 xlab=list(label=expression(paste(symbol("\341"),italic(l[vesicle]),symbol("\361"))),rot=30),
1227 ylab=list(label=expression(italic(l[monomer])),rot=-35),
1228 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1231 just=c("left","top"),
1232 gp=gpar(fontsize=48),
1233 y=unit(1,"npc")-unit(0.25,"lines"),
1234 x=unit(0,"npc")+unit(0.25,"lines"))
1237 \caption{Change in length backwards ($l_\mathrm{b}\times 10^{-6}$) (A) or activation energy
1238 in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the length of the
1239 monomer leaving ($l_\mathrm{monomer}$) and the average length of
1240 lipids in the vesicle ($\left<l_\mathrm{vesicle}\right>$).}
1241 \label{fig:lb_graph}
1245 \subsubsection{Complex Formation Backward}
1247 Complex formation ($CF1$) describes the interaction between \ac{CHOL} and
1248 \ac{PC} or \ac{SM}, where \ac{PC} or \ac{SM} protects the hydroxyl group of \ac{CHOL} from
1249 interactions with water %
1250 \citep{Huang1999:chol_solubility_model,
1251 Huang1999:cholesterol_solubility, McConnell2006, McConnell2003}%
1252 . \ac{PC} ($CF1=2$) can interact with two \ac{CHOL}, and \ac{SM} ($CF1=3$) with three
1253 \ac{CHOL} ($CF1=-1$). If the average of $CF1$ is positive (excess of \ac{SM} and
1254 \ac{PC} with regards to complex formation), components with negative $CF1$
1255 (\ac{CHOL}) will be retained. If average $CF1$ is negative, components with
1256 positive $CF1$ are retained. An equation which has this property is
1257 $CF1_\mathrm{b}=a^{\left<CF1_\mathrm{vesicle}\right>
1258 CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{vesicle}\right>
1259 CF1_\mathrm{monomer}\right|}$, where difference of the exponent is
1260 zero if the average $CF1$ and the $CF1$ of the component have the same
1261 sign, or double the product if the signs are different. Based on
1262 activation energy considerations, we took an ad hoc base for $a$ as
1267 CF1_\mathrm{b}=1.5^{\left<CF1_\mathrm{vesicle}\right> CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{vesicle}\right> CF1_\mathrm{monomer}\right|}
1268 \label{eq:complex_formation_backward}
1271 The most common $\left<CF1_\mathrm{vesicle}\right>$ is around $0.925$,
1272 which leads to a range of $\Delta \Delta G^\ddagger$ from
1273 $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*-1-abs(0.925*-1))))}
1274 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
1276 $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*2-abs(0.925*2))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
1277 for monomers with complex formation $2$ to
1278 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
1279 formation $0$. The full range of possible values for $CF1_\mathrm{b}$ are
1280 depicted in \cref{fig:cf1b_graph}.
1285 <<cf1b_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1286 trellis.device(new=F,color=TRUE)
1287 trellis.par.set(list(axis.line =list(col="transparent")))
1288 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1289 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1290 grid <- expand.grid(x=seq(-1,3,length.out=20),
1291 y=seq(-1,3,length.out=20))
1292 grid$z <- 1.5^(grid$x*grid$y-abs(grid$x*grid$y))
1293 print(wireframe(z~x*y,grid,cuts=50,
1296 scales=list(arrows=F,col=1,distance=0.75),
1297 xlab=list(label=expression(paste(symbol("\341"),italic(CF1[vesicle]),symbol("\361"))),rot=30),
1298 ylab=list(label=expression(italic(CF1[monomer])),rot=-35),
1299 zlab=list(label=expression(italic(CF1[b])),rot=93)),newpage=FALSE)
1302 just=c("left","top"),
1303 gp=gpar(fontsize=48),
1304 y=unit(1,"npc")-unit(0.25,"lines"),
1305 x=unit(0,"npc")+unit(0.25,"lines"))
1307 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1308 grid <- expand.grid(x=seq(-1,3,length.out=20),
1309 y=seq(-1,3,length.out=20))
1310 grid$z <- to.kcal(1.5^(grid$x*grid$y-abs(grid$x*grid$y)))
1311 print(wireframe(z~x*y,grid,cuts=50,
1314 scales=list(arrows=F,col=1,distance=0.75),
1315 xlab=list(label=expression(paste(symbol("\341"),italic(CF1[vesicle]),symbol("\361"))),rot=30),
1316 ylab=list(label=expression(italic(CF1[monomer])),rot=-35),
1317 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1320 just=c("left","top"),
1321 gp=gpar(fontsize=48),
1322 y=unit(1,"npc")-unit(0.25,"lines"),
1323 x=unit(0,"npc")+unit(0.25,"lines"))
1326 \caption{Change in complex formation 1 backwards ($CF1_\mathrm{b}$) (A) or
1327 activation energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus
1328 the complex formation of the monomer leaving
1329 ($CF1_\mathrm{monomer}$) and the average complex formation of the
1330 vesicle ($\left<CF1_\mathrm{vesicle}\right>$).}
1331 \label{fig:cf1b_graph}
1334 \section{Simulation Methodology}
1336 \subsection{Overall Architecture}
1338 The simulations are currently run using a program written in perl with
1339 various modules to handle the subsidiary parts. It produces output for
1340 each generation, including the step immediately preceding and
1341 immediately following a vesicle split (and optionally, each step) that
1342 is written to a state file which contains the state of the vesicle,
1343 environment, kinetic parameters, program invocation options, time, and
1344 various other parameters necessary to recreate the state vector at a
1345 given time. This output file is then read by a separate program in
1346 perl to produce different output which is generated out-of-band;
1347 output which includes graphs and statistical analysis is performed
1348 using R~\citep{R:maincite} (and various grid graphics modules) called
1349 from the perl program.
1351 The separation of simulation and output generation allows refining
1352 output, and simulations performed on different versions of the
1353 underlying code to be compared using the same analysis methods and
1354 code. It also allows later simulations to be restarted from a specific
1355 generation, utilizing the same environment.
1357 Random variables of different distributions are calculated using
1358 Math::Random (0.71), which is seeded for each run using entropy from
1359 the linux kernel's urandom device.
1361 Code is available upon request.
1363 \subsection{Environment Creation}
1365 \subsubsection{Components}
1366 The environment contains concentrations of components. In the current
1367 set of simulations, there are \Sexpr{1+4*5*7} different possible
1368 components, consisting of \ac{PC}, PE, \ac{PS}, \ac{SM}, and \ac{CHOL}, with all lipids
1369 except for \ac{CHOL} having 5 possible unsaturations ranging from 0 to 4,
1370 and 7 lengths from $12,14,...,22$ ($7\cdot
1371 5\cdot4+1=\Sexpr{1+4*5*7}$). In cases where the environment has less
1372 than the maximum number of components, the components are selected in
1373 order without replacement from a randomly shuffled deck of components
1374 (with the exception of \ac{CHOL}) represented once until the desired number
1375 of unique components is obtained. \ac{CHOL} is overrepresented
1376 $\Sexpr{5*7}$ times to be at the level of other lipid types, ensuring
1377 that the probability of \ac{CHOL} being absent in the environment is the
1378 same as the probability of any one of the other lipid types (\ac{PS}, PE,
1379 etc.) being absent. This reduces the number of simulations with a
1380 small number of components which are devoid of \ac{CHOL}.
1382 \subsubsection{Concentrations}
1383 Once the components of the environment have been selected, their
1384 concentrations are determined. In experiments where the environmental
1385 concentration is the same across all lipid components, the
1386 concentration is set to $10^{-10}\mathrm{M}$. In other cases, the
1387 environmental concentration is set to a random number from a gamma
1388 distribution with shape parameter 2 and an average of
1389 $10^{-10}\mathrm{M}$. The base concentration ($10^{-10}\mathrm{M}$)
1390 can also be altered at the initialization of the experiment to
1391 specific values for specific lipid types or components.
1393 The environment is a volume which is the maximum number of vesicles
1394 from a single simulation (4096) times the maximum size of the vesicle
1395 (usually 10000) divided by Avagadro's number divided by the total
1396 environmental lipid concentration, or usually
1397 \Sexpr{4096*10000/6.022E23/141E-10}~L.
1399 \subsection{Initial Vesicle Creation}
1401 Initial vesicles are seeded by selecting lipid molecules from the
1402 environment until the vesicle reaches a specific starting size. The
1403 vesicle starting size has gamma distribution with shape parameter 2
1404 and a mean of the per-simulation specified starting size, with a
1405 minimum of 5 lipid molecules. Lipid molecules are then selected to be
1406 added to the lipid membrane according to four different methods. In
1407 the constant method, molecules are added in direct proportion to their
1408 concentration in the environment. The uniform method adds molecules in
1409 proportion to their concentration in the environment scaled by a
1410 uniform random value, whereas the random method adds molecules in
1411 proportion to a uniform random value. The final method is a binomial
1412 method, which adjust the probability of adding a molecule of a
1413 specific component by the concentration of that component, and then
1414 adds the molecules one by one to the membrane. This final method is
1415 also used in the first three methods to add any missing molecules to
1416 the starting vesicle which were unallocated due to the requirement to
1417 add integer numbers of molecules. (For example, if a vesicle was to
1418 have 10 molecules split evenly between three lipid types, the 10th
1419 molecule would be assigned by randomly choosing between the three
1420 lipid types, weighted by their concentration in the environment.)
1422 \subsection{Simulation Step}
1424 Once the environment has been created and the initial vesicle has been
1425 formed, molecules join and leave the vesicle based on the kinetic
1426 parameters and state equation discussed above until the vesicle splits
1427 forming two progeny vesicles. The program then follows both vesicles
1428 and their progeny until the simulation reaches either the maximum
1429 number of vesicles (usually 4096), or the maximum simulation time
1432 \subsubsection{Calculation of Vesicle Properties}
1434 \label{sec:ves_prop}
1435 $S_\mathrm{vesicle}$ is the surface area of the vesicle, and is the sum of
1436 the surface area of all of the individual lipid components; each lipid
1437 type has a different surface area; we assume that the lipid packing is
1438 optimal, and there is no empty space.
1440 For the other vesicle properties (length, unsaturation, charge, and
1441 curvature), we calculate the mean, the standard deviation, the mean of
1442 the log, the mean of the absolute value of the log, the standard
1443 deviation of the log, and the standard deviation of the absolute value
1444 of the log. All cases, when we take the log, we take the log of the
1445 absolute value, and map $\log(0)$ to $0$. For the purpose of plotting
1446 vesicle properties, we only plot the mean of the property.
1448 \subsubsection{Joining and Leaving of Lipid Molecules}
1450 Determining the number of molecules to add to the lipid membrane
1451 ($n_i$) requires knowing $k_{\mathrm{f}i_\mathrm{adj}}$, the surface area of the
1452 vesicle $S_\mathrm{vesicle}$ (see \cref{sec:ves_prop}), the time interval
1453 $dt$ during which lipids are added, the base $k_{\mathrm{f}i}$, and the
1454 concentration of the monomer in the environment
1455 $\left[C_{i\mathrm{vesicle}}\right]$ (see \cref{eq:state}).
1456 $k_{\mathrm{f}i\mathrm{adj}}$ is calculated (see \cref{eq:kf_adj}) based on the
1457 vesicle properties and their hypothesized effect on the rate (in as
1458 many cases as possible, experimentally based)
1459 (see \cref{sec:kinetic_adjustments} for details). $dt$ can be varied
1460 (see \cref{sec:step_duration}), but for a given step is constant. This
1461 leads to the following:
1463 $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$
1465 In the cases where $n_i > 1$, the integer number of molecules is
1466 added. Fractional $n_i$ or the fractional remainder after the addition
1467 of the integer molecules are the probability of adding a specific
1468 molecule, and are compared to a uniformly distributed random value
1469 between 0 and 1. If the random value is less than or equal to the
1470 fractional part of $n_i$, an additional molecule is added.
1472 Molecules leaving the vesicle are handled in a similar manner, with
1474 $n_i = k_{\mathrm{b}i}k_{\mathrm{b}i_\mathrm{adj}}C_{i_\mathrm{vesicle}}\mathrm{N_A}dt$.
1476 Before addition or removal of molecules, the properties of the vesicle
1477 are calculated; this is done so that molecules can be considered to be
1478 added and removed at the same instant, even though additions are
1479 handled first programmatically. This also avoids cases where a removal
1480 would have resulted in a negative number of molecules for a particular
1483 \subsubsection{Step duration}
1484 \label{sec:step_duration}
1486 $dt$ is the time taken for each step of joining, leaving, and checking
1487 split. In the current implementation, it starts with a value of
1488 $10^{-6}$~$\mathrm{s}$ but this is modified in between each step if the
1489 number of molecules joining or leaving is too large or small. If more
1490 than half of the starting vesicle size molecules join or leave in a
1491 single step, $dt$ is reduced by half. If less than the starting
1492 vesicle size molecules join or leave in 100 steps, $dt$ is doubled.
1493 This is necessary to curtail run times and to automatically adjust to
1494 different experimental runs.
1496 \subsubsection{Vesicle splitting}
1498 If a vesicle has grown to a size which is more than double the
1499 starting vesicle size, the vesicle splits. Molecules are assigned to
1500 the progeny vesicles at random, with each progeny vesicle having an
1501 equal chance of getting a single molecule. The number of molecules to
1502 assign to the first vesicle has binomial distribution with a
1503 probability of an event in each trial of 0.5 and a number of trials
1504 equal to the number of molecules.
1508 The environment, initial vesicle, and the state of the vesicle
1509 immediately before and immediately after splitting are stored
1510 to produce later output.
1512 % \section{Analyzing output}
1514 % The analysis of output is handled by a separate perl program which
1515 % shares many common modules with the simulation program. Current output
1516 % includes simulation progress, summary tables, summary statistics, and
1519 % \subsection{PCA plots}
1521 % Two major groups of axes that contribute to vesicle variation between
1522 % subsequent generations are the components and properties of vesicles.
1523 % Each component in a vesicle is an axis on its own; it can be measured
1524 % either as an absolute number of molecules in each component, or the
1525 % fraction of molecules of that component over the total number of
1526 % molecules; the second approach is often more convenient, as it allows
1527 % vesicles with different number of molecules to be directly compared
1528 % (though it hides any effect of vesicle size).
1530 % In order to visualize the variation between and transition of
1531 % subsequent generations of vesicles from their initial state in the
1532 % simulation, to their final state at the termination of the simulation,
1533 % we plot the projection of the generations onto the two principle PCA
1534 % axes (and alternatively, any pairing of the axes).
1536 % \subsection{Carpet plots}
1538 % Carpet plots show the distance/similarity between the composition or
1539 % properties of all generations in a single run. The difference between
1540 % each group of vesicle can be calculated using Euclidean distance
1541 % (properties and compositions) or H similarity (composition only). We
1542 % must use Euclidean distance for properties because the H distance
1543 % metric is invalid when comparing negative vector elements to positive
1546 % In addition to showing distance or similarity, carpet plots also
1547 % depict propersomes and composomes as square boxes on the diagonals and
1548 % propertypes and compotypes as rectangles off the diagonals, each
1549 % propertype or compotype with a distinct color.
1551 % \subsection{Propersomes, propertypes, composomes and compotypes}
1553 % A generation is considered to be a propersome if it is less than $0.6$
1554 % units (by Euclidean distance of normalized properties) away from the
1555 % generation immediately following and preceding. Likewise, a generation
1556 % is in a composome if its H similarity is more than $0.9$ (by the
1557 % normalized H metric) from the preceding and following generations.
1558 % Propersomes and composomes are then assigned to propertypes and
1559 % compotypes using Paritioning Around Medioids
1560 % (PAM). All values of $k$ between 2 and 15
1561 % (or the number of propersomes and composomes, if that is less than 15)
1562 % are tried, and the clustering with the smallest
1563 % silhouette~\citep{Rousseeuw1987:silhouettes} is chosen as the ideal
1564 % clustering~\citep{Shenhav2005:pgard}.
1567 %\bibliographystyle{unsrtnat}
1568 %\bibliography{references.bib}