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 = -1, digits = 2)
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 (\cref{tab:kinetic_parameters_lipid_types}).
158 \begin{tabular}{c c c c c c c c}
160 Type & $k_\mathrm{f}$ $\left(\frac{\mathrm{m}}{\mathrm{s}}\right)$
161 & $k'_\mathrm{f}$ $\left(\frac{1}{\mathrm{M} \mathrm{s}}\right)$
162 & $k_\mathrm{b}$ $\left(\mathrm{s}^{-1}\right)$
163 & Area $\left({Å}^2\right)$ & Charge & $\mathrm{CF}1$ & Curvature \\
165 PC & $\Sexpr{kf[1]}$ & $3.7 \times 10^6$ & $2 \times 10^{-5}$ & 63 & 0 & 2 & 0.8 \\
166 PS & $\Sexpr{kf[2]}$ & $3.7 \times 10^6$ & $1.25\times 10^{-5}$ & 54 & -1 & 0 & 1 \\
167 CHOL & $\Sexpr{kf[3]}$ & $5.1 \times 10^7$ & $2.8 \times 10^{-4}$ & 38 & 0 & -1 & 1.21 \\
168 SM & $\Sexpr{kf[4]}$ & $3.7 \times 10^6$ & $3.1 \times 10^{-3}$ & 61 & 0 & 3 & 0.8 \\
169 PE & $\Sexpr{kf[5]}$ & $2.3 \times 10^6$ & $1 \times 10^{-5}$ & 55 & 0 & 0 & 1.33 \\
172 \caption{Kinetic parameters and molecular properties of lipid types}
173 \label{tab:kinetic_parameters_lipid_types}
176 %%% \DLA{I think we may just reduce these three sections; area, $k_\mathrm{f}$
177 %%% and $k_\mathrm{b}$ to \cref{tab:kinetic_parameters_lipid_types} with
180 %%% \RZ{Yes, but we also have to have then as comments the numbers that
181 %%% are not supported by refs}
183 \subsubsection{$k_\mathrm{f}$ for lipid types}
184 $k_{\mathrm{f}_\mathrm{PC}}$ was measured by
185 \citet{Nichols1985:phospholipid_monomer_vesicle_thermodynamics} and
186 was found to be $3.7\times 10^6$~$\frac{1}{\mathrm{M}\mathrm{s}}$ by the
187 partitioning of \ac{NBDPC} between \ac{DOPC} vesicles and water. As
188 similar references are not available for \ac{SM} or \ac{PS}, we assume that they have
189 the same $k_\mathrm{f}$. For \ac{CHOL}, no direct measurement of $k_\mathrm{f}$ is available,
190 however, \citet{Estronca2007:dhe_kinetics} measured the transfer of
191 \ac{DHE} from \ac{BSA} to \acp{LUV}, and found a $k_\mathrm{f}$ of $5.1\times 10^7$~%
192 $\frac{1}{\mathrm{M} \mathrm{s}}$. We assume that this value is close
193 to that of \ac{CHOL}, and use it for $k_{\mathrm{f}_\mathrm{\ac{CHOL}}}$. In the case of
194 \ac{PE}, \citet{Abreu2004:kinetics_ld_lo} measured the association of
195 \ac{NBDDMPE} with \ac{POPC} \acp{LUV} and found a value for $k_\mathrm{f}$ of $2.3 \times 10^{6}$~%
196 $\frac{1}{\mathrm{M} \mathrm{s}}$. These three authors used a slightly
197 different kinetic formalism ($\frac{d\left[A_v\right]}{dt} =
198 k'_\mathrm{f}[A_m][L_v] - k_\mathrm{b}[A_v]$), so we correct their values of $k_\mathrm{f}$ by
199 multiplying by the surface area of a mole of lipids.
201 \subsubsection{$k_\mathrm{b}$ for lipid types}
203 \citet{Wimley1990:dmpc_exchange} measured the half time for the
204 exchange of [$^3$H]\ac{DMPC} between \acp{LUV} at 30\textdegree C, and found it
205 to be 9.6 hr. As this is a first order reaction, and the primary
206 limiting step in exchange is lipid desorption, $k_\mathrm{b}$ for \ac{DMPC} is
207 $k_{\mathrm{b}_\mathrm{PC}}=\frac{\log 2}{9.6 \times 60 \times 60} \approx
208 \Sexpr{log(2)/(9.6*60*60)}
209 \mathrm{s}^{-1}$. We assume that $k_\mathrm{b}$ for \ac{SM} is the same as for \ac{PC}.
210 To estimate the $k_\mathrm{b}$ of \ac{PE} and \ac{PS}, we used the data from
211 \citet{Nichols1982:ret_amphiphile_transfer} who measured the rate of
212 exchange of the fluorescent label \ac{C6NBD} attached to different
213 lipid species. Although the values of $k_\mathrm{b}$ are different for the labeled
214 and unlabeled lipids, we assume that the ratios of the kinetics
215 constants for the lipid types are the same. Furthermore we assume that
216 PG behaves similarly to \ac{PS}. Thus, we can determine the $k_\mathrm{b}$ of \ac{PE} and
217 \ac{PS} from the already known $k_\mathrm{b}$ of \ac{PC}. For \ac{C6NBD} labeled \ac{PC},
218 \citet{Nichols1982:ret_amphiphile_transfer} obtained a $k_\mathrm{b}$ of
219 $0.89$~$\mathrm{min}^{-1}$, \ac{PE} of $0.45$~$\mathrm{min}^{-1}$ and PG of
220 $0.55$~$\mathrm{min}^{-1}$. Assuming the ratios of the rate of
221 exchange are the same for unlabeled lipids and labeled lipids, we can
222 determine the $k_\mathrm{b}$ of \ac{PE} and \ac{PS} from the already known $k_\mathrm{b}$ of
223 \ac{PC}~\citep{Wimley1990:dmpc_exchange}. Calculating the ratio, this leads
224 us to $k_{\mathrm{b}_\mathrm{PE}} =
225 \frac{k_{\mathrm{b}_\mathrm{PC}}\times\mathrm{PE}}{\mathrm{PC}} \approx
226 \frac{2\times 10^{-5}\,\mathrm{s}^{-1} \times
227 0.45\,\mathrm{min}^{-1}}{0.89\,\mathrm{min}^{-1}} \approx
228 \Sexpr{log(2)/(9.6*60*60)*0.45/0.89}$~$\mathrm{s}^{-1}$
229 and likewise, $k_{\mathrm{b}_\mathrm{PS}}\approx
230 \Sexpr{log(2)/(9.6*60*60)*0.55/0.89}$~$\mathrm{s}^{-1}$.
231 The $k_\mathrm{b}$ of \ac{SM} was determined using the work of
232 \citet{Bai1997:lipid_movementbodipy}, who measured spontaneous
233 transfer of C$_5$-DMB-SM and C$_5$-DMB-PC from donor and acceptor
234 vesicles, finding $3.4\times10^{-2}$~$\mathrm{s}^{-1}$ and
235 $2.2\times10^{-3}$~$\mathrm{s}^{-1}$ respectively; using the ratio of
236 $k_\mathrm{b}$ of C$_5$-DMB-SM to the $k_\mathrm{b}$ of C$_5$-DMB-PC times the $k_\mathrm{b}$ of
237 \ac{PC} ($\frac{3.4 \times 10^{-2}\mathrm{s}^{-1}}{2.2 \times
238 10^{-3}\mathrm{s}^{-1}}\approx
239 \Sexpr{log(2)/(9.6*60*60)}\mathrm{s}^{-1}$),
240 we obtain $k_{\mathrm{b}_\mathrm{SM}} \approx
241 \Sexpr{log(2)/(9.6*60*60)*3.4e-2/2.2e-3}$.
242 In the case of \ac{CHOL}, \citet{Jones1990:lipid_transfer} measured the
243 $t_{1/2}$ of [$^3$H] transfer from POPC vesicles and found it to be 41
244 minutes, leading to a $k_{\mathrm{b}_\mathrm{CHOL}} = \frac{\log 2}{41\times
246 \Sexpr{log(2)/(41*60)}$~$\mathrm{s}^{-1}$.
248 \subsubsection{Headgroup Surface Area for lipid types}
250 % Smondyrev, A. M., and M. L. Berkowitz. 1999b. United atom force
251 % field for phospholipid membranes: constant pressure molecular
252 % dynamics simulation of dipalmitoylphosphatidicholine/water system.
253 % J. Comput. Chem. 50:531–545.
255 Different lipids have different headgroup surface areas, which contributes to
256 $\left[S_\mathrm{vesicle}\right]$. \citet{Smaby1997:pc_area_with_chol}
257 measured the surface area of POPC with a Langmuir film balance, and
258 found it to be 63~Å$^2$ at $30$~$\frac{\mathrm{mN}}{\mathrm{m}}$.
259 Molecular dynamic simulations found an area of 54 Å$^2$ for
260 DPPS\citep{Cascales1996:mds_dpps_area,Pandit2002:mds_dpps}, which is
261 in agreement with the experimental value of 56~Å$^2$ found using a
262 Langmuir balance by \citet{Demel1987:ps_area}.
263 \citet{Shaikh2002:pe_phase_sm_area} measured the area of \ac{SM} using a
264 Langmuir film balance, and found it to be 61~Å$^2$. Using $^2$H NMR,
265 \citet{Thurmond1991:area_of_pc_pe_2hnmr} found the area of
266 DPPE-d$_{62}$ to be 55.4 Å$^2$. \citet{Robinson1995:mds_chol_area}
267 found an area for \ac{CHOL} of 38~Å$^2$ using molecular dynamic
270 % robinson's chol area is kind of crappy; they did it using MDS, but
271 % vaguely; they don't indicate what it was at the end, and they only
272 % had 2 molecules of \ac{CHOL}.
276 % \ac{PC}: 63 A (Smaby97rd) (Pandit (?)
277 % \ac{PS}: 54 A (Pandit02LIB) (Cascales 1996; J. Chem. Phys 104:2713, Mukhopadhyay2004)
278 % \ac{CHOL}: 38 A (Robinson 1995) (Lundberg 1982)
279 % \ac{SM}: 51 A (Shaikh2002, but 61(?))
280 % PE: 55 A (Thurmond, 1991) (Pandit2002)
282 % Compare to results by Petrache2004 with DOPC of 72.5, DOPS of 65.3.
284 % The cross‑sectional area of cholesterol is ~37 ֵ2 (Lundberg, B. 1982. A
285 % surface film study of the lateral packing of phosphatidylcholine and
286 % cholesterol. Chem. Phys. Lipids. 31:23-32).
288 % From Pandit02LIB: the area per headgroup for PE is ~10‑20\% (Thurmond
289 % et al., 1991) smaller than the area per PC molecule. For DLPE, it is
290 % 50.6 ֵ2 (McIntosh and Simon, 1986; Damodaran et al., 1992). The
291 % estimated area per molecule for DPPE is ~55.4 ֵ2 (Thurmond et al.,
294 % From Pandit02LIB: (from their simulations) average area per headgroup
295 % for DPPS molecules is 53.75 ± 0.10 ֵ2. The values inferred from the
296 % experiments are in the region between 45 and 55 ֵ2 (Cevc et al., 1981.
297 % Titration of the phase transition of phosphatidylserine bilayer
298 % membranes. Effect of pH, surface electrostatics, ion binding, and
299 % head‑group hydration. Biochemistry. 20:4955‑4965; | Demel et al.,
300 % 1987. Monolayer characteristics and thermal behavior of natural and
301 % synthetic phosphatidylserines. Biochemistry. 26:8659‑8665). Note that
302 % the average area per headgroup for DPPC bilayer obtained from
303 % simulations and measurements is ~62 ֵ A2. One would expect that the
304 % area per headgroup in case of DPPS molecules will be larger than the
305 % area per DPPC molecule, because the headgroups of DPPS are charged and
306 % therefore repel each other. Contrary to this expectation, the area per
307 % headgroup in DPPS is ~13\% smaller than that of DPPC. López Cascales
308 % et al. (1996. Molecular dynamics simulation of a charged biological
309 % membrane. J. Chem. Phys. 104:2713‑2720.) proposed that this reduction
310 % in the area per headgroup is due to a strong intermolecular
311 % coordination between DPPS molecules. MD gave area per POPS of 55 A2 at
312 % 300K (Mukhopadhyay04LIB).
314 % However, by 2H NMR and X-ray, Petrache04LIB determined the area of
315 % DOPS at 30 C to be 65.3 A2, considerably smaller than DOPC, determined
316 % by the same group in another paper to be 72.5 ֵ A2
317 % (Tristram‑Nagle00LIB).
319 % From Shaikh02: At 35 mN/m of surface pressure, the areas/molecule for
320 % POPE, SM, and CHOL were found to be 65.2, 61.1, and 38.8 Å2.
322 % In Lonchin99: POPC area 72 A2 ((a) Cornell, B. A.; Middlehurst, J.;
323 % Separovic, F. Biochim. Biophys. Acta1980, 598, 405−410), oleic acid
324 % area 32 A2 ((a) Small, D. M. In Handbook of Lipid Research; Hanahan,
325 % D. J., Ed.; Plenum Press: New York, 1986; Vol. 4, Chapter 9.
327 % (b) Jain, M. K. Introduction to Biological Membranes, 2nd ed.; John
328 % Wiley \& Sons: New York, 1988; Chapter 3). Cistola, D. P.; Atkinson,
329 % D.; Hamilton, J. A.; Small, D. M. Biochemistry1986, 25, 2804−2812.
331 % From Marsh96: The values obtained for the headgroup area per molecule,
332 % fitted as free parameters in Fig. 6, are APE = 58 A2 and APC = 78 A2,
333 % for DOPE and DOPC, respectively. These values are within the range
334 % obtained for the corresponding areas per molecule of
335 % phosphatidylethanolamines and phosphatidylcholines in the fluid
336 % lamellar phase (see, e.g., Marsh, 1990), and correspond to values for
337 % the lipid packing parameter of V/Al = 1.08 and 1.34 for DOPC and DOPE,
338 % respectively. The former value is consistent with DOPC alone forming
339 % hydrated lamellar structures rather than nonlamellar structures. The
340 % latter value lies at the lower end of the packing parameters obtained
341 % for phosphatidylethanolamines.
343 % From Kumar91: The head group area chosen was 71.7 A2 for PC and
344 % lysoPC, 42 A2 for PE, and 19 A2 for Chol (5, 14). The value selected
345 % for PC is in agreement with the 71 A2 value determined for PCs with
346 % saturated acyl chains (14). The hydrocarbon volume and the length
347 % chosen for Chol are 400 A3 and 17.25 A, respectively (5).
351 % From Sampaio05: Besides this work and our own earlier report on the
352 % association of NBD-DMPE with lipid bilayers (Abreu et al., 2004), we
353 % are aware of only one other report in the literature (Nichols, 1985)
354 % in which all the kinetic constants of lipid-derived amphiphile
355 % association with lipid bilayer membranes were experimentally measured.
356 % {indeed; everything is k- !!!; rz}
358 % From McLean84LIB: Although it is difficult to measure cmc values for
359 % the sparingly soluble lipids used in this study, estimates for
360 % lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds,
361 % 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution;
362 % Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6
363 % X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X
364 % 10-8 M for DMPC was estimated from the linear relationship between ln
365 % cmc and the number of carbons in the PC acyl chain by using data for n
366 % = 7, 8, 10, and 16 [summarized in Tanford (1980)].
368 % From Toyota08: Recently, several research groups have reported
369 % self-reproducing systems of giant vesicles that undergo a series of
370 % sequential transformations: autonomous growth, self-division, and
371 % chemical reactions to produce membrane constituents within the giant
374 % Vesicle sizes of 25 nm for SUV and 150 nm for LUV were mentioned by
377 % From Lund-Katz88: Charged and neutral small unilamellar vesicles
378 % composed of either saturated PC, unsaturated PC, or SM had similar
379 % size distributions with diameters of 23 \& 2 nm.
381 % From Sampaio05LIB: The exchange of lipids and lipid derivatives
382 % between lipid bilayer vesicles has been studied for at least the last
383 % 30 years. Most of this work has examined the exchange of amphiphilic
384 % molecules between a donor and an acceptor population. The measured
385 % efflux rates were shown in almost all cases, not surprisingly, to be
386 % first order processes. In all of this work, the insertion rate has
387 % been assumed to be much faster than the efflux rate. Having measured
388 % both the insertion and desorption rate constants for amphiphile
389 % association with membranes, our results show that this assumption is
390 % valid. In several cases reported in the literature, the insertion rate
391 % constant was assumed, although never demonstrated, to be a
392 % diffusion-controlled process.
394 % (for methods? From McLean84LIB: The activation free energies and free
395 % energies of transfer from self-micelles to water increase by 2.2 and
396 % 2.1 kJ mol-' per methylene group, respectively. {see if we can use it
397 % to justify arranging our changed in activating energy around 1
400 % Jones90 give diameter of LUV as 100 nm, and of SUV as 20 nm; that
401 % would give the number of molecules per outer leaflet of a vesicle as
404 % Form Simard08: Parallel studies with SUV and LUV revealed that
405 % although membrane curvature does have a small effect on the absolute
406 % rates of FA transfer between vesicles, the ΔG of membrane desorption
407 % is unchanged, suggesting that the physical chemical properties which
408 % govern FA desorption are dependent on the dissociating molecule rather
409 % than on membrane curvature. However, disagreements on this fundamental
410 % issue continue (57, 61, 63, 64)
412 % (methods regarding the curvature effect: Kleinfeld93 showed that the
413 % transfer parameters of long-chain FFA between the lipid vesicles
414 % depend on vesicle curvature and composition. Transfer of stearic acid
415 % is much slower from LUV as compared to SUV).
417 % From McLean84: In a well-defined experimental system consisting of
418 % unilamellar lipid vesicles, in the absence of protein, the
419 % rate-limiting step for the overall exchange process is desorption
420 % (McLean \& Phillips, 1981). {thus I can take exchange data for the
421 % estimation of k- rz; 8/11/08}.
423 \subsubsection{Complex Formation 1}
425 \citet{Thomas1988:chol_transfer} found that \ac{SM} significantly decreases
426 the rate of cholesterol transfer from \ac{PC}, while \ac{PE} and \ac{PS} have no
427 effect at physiologically significant levels. We attribute this to the
428 formation of complexes between \ac{SM} and \ac{CHOL}, which retards the rate of
429 efflux of \ac{CHOL} from the membrane. Similar complexes exist between \ac{PC}
430 and \ac{CHOL}, where it was shown that two \ac{CHOL} molecules complex with one
431 \ac{PC}~\citep{Huang1999:chol_solubility_model,
432 Huang1999:cholesterol_solubility,McConnell2006,McConnell2003}. It
433 has also been shown that \ac{SM} binds more \ac{CHOL} molecules than
434 \ac{PC}~\citep{Katz1988:pl_packing_chol}. We assume that three \ac{CHOL} complex
435 with one \ac{SM}, leading to $\mathrm{CF}1$ values of $3$, $2$, and $-1$
436 for \ac{SM}, \ac{PC}, and \ac{CHOL}, respectively.
438 \subsubsection{Curvature}
440 We used the formulation for curvature of
441 \citet{Israelachvili1975:amphiphile_self_assembly}, or $S=\frac{v}{l_c
442 a}$, where $S$ is the curvature, which ranges from $0$ to $\infty$,
443 $l_c$ is the critical length of the acyl chain, $v$ is the volume of
444 the lipid, and $a$ is the area of the head group.
445 \citet{Kumar1991:lipid_packing} found a value for $S$ of 1.33 for PE,
446 1.21 for \ac{CHOL}, and 0.8 for diC$_{16}$ \ac{PC}. We assume that \ac{PS} has neutral
447 curvature (value of 1), and that \ac{SM} has the same curvature as \ac{PC}
448 (0.8). It should also be noted that in reality the curvature parameter
449 changes with length, but at longer chain lengths, is effectively
450 constant; in the current model, curvature is held constant for each
454 % 1. Israelachvili, J. N., Mitchell, D. J. & Ninham, B. W. (1976) J.
455 % Chem. Soc. Faraday Trans. 2 72, 1525-1568.
456 % 2. Israelachvili, J. N., Marcelja, S. & Horn, R. G. (1980) Q. Rev.
457 % Biophys. 13, 121-200.
458 % 3. Israelachvili, J. N. (1985) Intermolecular and Surface Forces
459 % (Academic, New York), pp. 229-271.
461 % \DLA{This section needs to be added still.}
463 % PE = 1.33 (Kumar91)
464 % CHOL = 1.21 (Kumar91)
466 % SM=0.8 (assumed by rz same as PC)
467 % PS=1 (no refs so far; should be close to unity; rz)
469 % From Epand94LIB: a considerable amount of heat is evolved per mole of
470 % bilayer stabilizer when such molecules are inserted into membranes
471 % which are under a high curvature strain. If this energy were coupled
472 % to a membrane event, such as the conformational change in a protein,
473 % it could be the driving force responsible for such processes.
475 % From Epand94LIB: considerable energy may be released upon the
476 % incorporation of certain molecules into membranes which have a low
477 % radius of spontaneous curvature. {discuss with DL in terms of
478 % catalytic activity; rz; 8/24/2010}.
480 % Kleinfeld93 showed that the transfer parameters of long-chain FFA
481 % between the lipid vesicles depend on vesicle curvature and
482 % composition. Transfer of stearic acid is much slower from LUV as
485 % Substitution of DPPC SUV for POPC SUV as the donor matrix revealed
486 % some differences in the energetics of transfer that are associated
487 % with the physical state of the donor particle. Whereas the free
488 % energies of transfer of oleic acid from DPPC SUV and POPC SUV are
489 % similar, the free energy of transfer of stearic acid from DPPC SUV is
490 % much higher than that observed from POPC.
492 % From Kumar91: S is given by S = V/al, where V is the hydrocarbon
493 % volume, a is the area of head group, and l is the critical length of
494 % the hydrocarbon chain (1-3). a, V, and I are all estimable or
495 % measurable (4), and the value of S can be calculated. The value of S
496 % determines the aggregate formed by lipids or any amphiphiles upon
497 % hydration. It has been shown that lipids aggregate to form spherical
498 % micelles (S < 1/3), nonspherical (cylindrical) micelles (1/3 < S <
499 % 1/2), bilayers (1/2 < S < 1), and reverse micelles or hexagonal (HII,)
500 % phases (S > 1). However, theoreticians caution that the above
501 % predicted limits set on the values of S are relatively insensitive to
502 % the exact values of V and a but are strongly dependent upon the choice
505 % From Kumar91: Cylindrical and wedge-shaped molecules have been shown
506 % to be essential for spontaneous vesiculation (6) as well as bilayer
507 % stabilization (7, 8).
509 % From Kumar91: It has also been shown experimentally that molecules
510 % having complementary molecular shapes can form bilayer structures.
511 % Inverted cone (wedge)- shaped lysoPC and cone-shaped Chol or
512 % unsaturated PE can interact stoichiometrically to form bilayer
513 % structures (9-12). {quote to justify our formula for complementary
516 \section{Kinetic adjustments}
517 \label{sec:kinetic_adjustments}
519 In the mass action equation used in the formalism, both the forward
520 and backwards kinetic constants ($k_\mathrm{f}$ and $k_\mathrm{b}$) are adjusted (by
521 $k_{\mathrm{f}i\mathrm{adj}}$ and $k_{\mathrm{b}i\mathrm{adj}}$) to reflect the
522 properties of the vesicle and the properties of the monomer entering
523 or exiting the vesicle.
525 \subsection{Forward Kinetic Adjustments}
527 The forward rate constant adjustment, $k_{\mathrm{f}i\mathrm{adj}}$ takes into
528 account unsaturation ($un_\mathrm{f}$), charge ($ch_\mathrm{f}$), curvature ($cu_\mathrm{f}$),
529 length ($l_\mathrm{f}$), and complex formation ($CF1_\mathrm{f}$), each of which is
530 modified depending on the specific component and the vesicle into
531 which the component is entering; the final $k_{\mathrm{f}i\mathrm{adj}}$ is the
532 product of the adjustments made separately for each property.
535 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}
539 \subsubsection{Unsaturation Forward}
541 In order for a lipid to be inserted into a membrane, a void has to be
542 formed for it to fill. Voids can be generated by the combination of
543 unsaturated and saturated lipids forming heterogeneous domains. Void
544 formation is increased when the unsaturation of lipids in the vesicle
545 is widely distributed; in other words, the insertion of lipids into
546 the membrane is greater when the standard deviation of the
547 unsaturation is larger. %
548 %%% \RZ{May need to site (at least for us) works showing
549 %%% mismatch-dependent ``defects''}. %
550 Assuming that an increase in width of the distribution linearly
551 decreases the free energy of activation, the $un_\mathrm{f}$ parameter must
552 follow $a^{\mathrm{stdev}\left(un_\mathrm{vesicle}\right)}$ where $a > 1$,
553 so a convenient starting base for $a$ which leads to a reasonable
554 change in Eyring activation energy is $2$:
557 un_\mathrm{f} = 2^{\mathrm{stdev}\left(un_\mathrm{vesicle}\right)}
558 \label{eq:unsaturation_forward}
561 The mean $\mathrm{stdev}\left(un_\mathrm{vesicle}\right)$ in our
562 simulations is around $1.5$, which leads to a $\Delta \Delta
563 G^\ddagger$ of $\Sexpr{to.latex(format(digits=3,to.kcal(2^1.5)))}
564 \frac{\mathrm{kcal}}{\mathrm{mol}}$, and a total range of possible
565 values depicted in \cref{fig:unf_graph}.
567 % \RZ{Explain here, or even earlier that the formulas were ad hoc
568 % adjusted to correspond to ``reasonable'' changes in the Eyring
569 % activation energies.}
571 % It is not clear that the unsaturation of the inserted monomer will
572 % affect the rate of the insertion positively or negatively, so we do
573 % not include a term for it in this formalism.
577 <<unf_graph,echo=FALSE,results='hide',fig.width=10,fig.height=5,cache=FALSE>>=
578 layout(matrix(1:2,nrow=1,ncol=2))
579 curve(2^x,from=0,to=sd(c(0,4)),
580 xlab=expression(paste(stdev,group("(",italic(un[vesicle]),")"))),
581 ylab=expression(italic(un[f])))
582 vps <- baseViewports()
583 pushViewport(vps$figure)
585 just=c("left","top"),
586 gp=gpar(fontsize=48),
587 y=unit(1,"npc")-unit(0.25,"lines"),
588 x=unit(0,"npc")+unit(0.25,"lines"))
590 curve(to.kcal(2^x),from=0,to=sd(c(0,4)),
591 xlab=expression(paste(stdev,group("(",italic(un[vesicle]),")"))),
592 ylab="Activation Energy (kcal/mol)")
593 vps <- baseViewports()
594 pushViewport(vps$figure)
596 just=c("left","top"),
597 gp=gpar(fontsize=48),
598 y=unit(1,"npc")-unit(0.25,"lines"),
599 x=unit(0,"npc")+unit(0.25,"lines"))
602 \caption{Change in unsaturation forward ($un_\mathrm{f}$) (A) or activation
603 energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the
604 standard deviation of the unsaturation of the vesicle
605 ($\mathrm{stdev}\left(un_\mathrm{vesicle}\right)$).}
606 \label{fig:unf_graph}
609 \subsubsection{Charge Forward}
611 A charged lipid such as \ac{PS} approaching a vesicle with an average
612 charge of the same sign will experience repulsion, whereas those with
613 different signs will experience attraction, the degree of which is
614 dependent upon the charge of the monomer and the average charge of the
615 vesicle. If either the vesicle or the monomer has no charge, there
616 should be no effect of charge upon the rate. This leads us to the
617 following equation, $a^{-\left<ch_v\right> ch_m}$ (which is similar to
618 the numerator of Coulomb's law), where $\left<ch_v\right>$ is the
619 average charge of the vesicle, and $ch_m$ is the charge of the
620 monomer. If either $\left<ch_v\right>$ or $ch_m$ is 0, the adjustment
621 parameter is 1 (no change), whereas it decreases if both are positive
622 or negative, as the product of two real numbers with the same sign is
623 always positive. The base for $a$ ($60$) was chosen ad hoc to
624 correspond to a maximum of about
625 $0.5\frac{\mathrm{kcal}}{\mathrm{mol}}$ change in activation energy
626 for the common case, resulting in the following equation:
630 ch_\mathrm{f} = 60^{-\left<{ch}_v\right> {ch}_m}
631 \label{eq:charge_forward}
634 The most common $\left<{ch}_v\right>$ is around $-0.165$, which leads
635 to a range of $\Delta \Delta G^\ddagger$ from
636 $\Sexpr{format(digits=3,to.kcal(60^(-.165*-1)))}
637 \frac{\mathrm{kcal}}{\mathrm{mol}}$ to
638 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$, and the total range of possible
639 values seen in \cref{fig:chf_graph}.
643 <<chf_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
644 trellis.device(new=F,color=TRUE)
645 trellis.par.set(list(axis.line =list(col="transparent")))
646 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
647 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
648 x <- seq(-1,0,length.out=20)
649 y <- seq(-1,0,length.out=20)
650 grid <- expand.grid(x=x,y=y)
651 grid$z <- as.vector(60^(-outer(x,y)))
652 print(wireframe(z~x*y,grid,cuts=50,
655 scales=list(arrows=F,col=1,distance=0.75),
656 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
657 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
658 zlab=list(label=expression(italic(ch[f])),rot=93)),newpage=FALSE)
661 just=c("left","top"),
662 gp=gpar(fontsize=48),
663 y=unit(1,"npc")-unit(0.25,"lines"),
664 x=unit(0,"npc")+unit(0.25,"lines"))
666 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
667 x <- seq(-1,0,length.out=20)
668 y <- seq(-1,0,length.out=20)
669 grid <- expand.grid(x=x,y=y)
670 grid$z <- as.vector(to.kcal(60^(-outer(x,y))))
671 print(wireframe(z~x*y,grid,cuts=50,
674 scales=list(arrows=F,col=1,distance=0.75),
675 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
676 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
677 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
680 just=c("left","top"),
681 gp=gpar(fontsize=48),
682 y=unit(1,"npc")-unit(0.25,"lines"),
683 x=unit(0,"npc")+unit(0.25,"lines"))
686 \caption{Change in charge forward ($ch_\mathrm{f}$) (A) or activation energy in
687 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the charge of the
688 monomer entering ($ch_\mathrm{monomer}$) and the average charge of
689 the vesicle ($\left<ch_\mathrm{vesicle}\right>$).}
690 \label{fig:chf_graph}
693 \subsubsection{Curvature Forward}
695 Curvature is a measure of the intrinsic propensity of specific lipids
696 to form micelles (positive curvature), inverted micelles (negative
697 curvature), or planar sheets (neutral
698 curvature)~\citep{Israelachvili1975:amphiphile_self_assembly}. %
699 In this formalism, curvature is measured as the ratio of the volume of
700 the lipid to the area of the head times the length of the lipid
701 ($S=\frac{v}{l_c a}$), so negative curvature is bounded by $(0,1)$,
702 neutral curvature is 1, and positive curvature is bounded by
703 $(1,\infty)$. The curvature can be transformed using $\log$, which has
704 the property of making the range of positive and negative curvature
705 equal, and distributed about 0.
707 As in the case of unsaturation, void formation is increased by the
708 presence of lipids with mismatched curvature. Thus, a larger
709 distribution of curvature in the vesicle increases the rate of lipid
710 insertion into the vesicle. However, a component with curvature
711 $cu^{-1}$ will cancel out a component with curvature $cu$, so we have
712 to log transform (turning these into $-\log cu$ and $\log cu$), then
713 take the absolute value ($\log cu$ and $\log cu$), and finally measure
714 the width of the distribution, which in the case of exactly mismatched
715 curvatures is $0$. Thus, by using the log transform to make the range
716 of the lipid curvature equal between positive and negative, and taking
717 the average to cancel out exactly mismatched curvatures, we come to an
718 equation with the shape $a^{\left<\log cu_\mathrm{vesicle}\right>}$.
719 An ad hoc base for $a$ which for the common case is $10$, yielding:
723 % cu_\mathrm{f} = 10^{\mathrm{stdev}\left|\log cu_\mathrm{vesicle}\right|}
724 cu_\mathrm{f} = 10^{\left|\left<\log cu_\mathrm{monomer} \right>\right|\mathrm{stdev} \left|\log cu_\mathrm{vesicle}\right|}
725 \label{eq:curvature_forward}
728 The most common $\left|\left<\log {cu}_v\right>\right|$ is around
729 $0.013$, which with the most common $\mathrm{stdev} \log
730 cu_\mathrm{vesicle}$ of $0.213$ leads to a $\Delta \Delta G^\ddagger$
731 of $\Sexpr{format(digits=3,to.kcal(10^(0.13*0.213)))}
732 \frac{\mathrm{kcal}}{\mathrm{mol}}$. This is a consequence of the
733 relatively matched curvatures in our environment. The full range of
734 $cu_\mathrm{f}$ values possible are shown in \cref{fig:cuf_graph}.
736 % 1.5 to 0.75 3 to 0.33
738 <<cuf_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
739 trellis.device(new=F,color=TRUE)
740 trellis.par.set(list(axis.line =list(col="transparent")))
741 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
742 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
743 grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
744 sd(abs(log(c(1,0.33)))),
745 sd(abs(log(c(0.33,3)))))),length.out=20),
746 y=seq(0,max(c(mean(log(c(1,3)),
747 mean(log(c(1,0.33))),
748 mean(log(c(0.33,3)))))),length.out=20))
749 grid$z <- 10^(grid$x*grid$y)
750 print(wireframe(z~x*y,grid,cuts=50,
753 scales=list(arrows=F,col=1,distance=0.75),
754 xlab=list(label=expression(paste(stdev,group("(",italic(log(cu[vesicle])),")"))),rot=30),
755 ylab=list(label=expression(paste(symbol("\341"),italic(log(cu[vesicle])),symbol("\361"))),rot=-35),
756 zlab=list(label=expression(italic(cu[f])),rot=93)),newpage=FALSE)
759 just=c("left","top"),
760 gp=gpar(fontsize=48),
761 y=unit(1,"npc")-unit(0.25,"lines"),
762 x=unit(0,"npc")+unit(0.25,"lines"))
764 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
765 grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))),
766 sd(abs(log(c(1,0.33)))),sd(abs(log(c(0.33,3)))))),length.out=20),
767 y=seq(0,max(c(mean(log(c(1,3)),
768 mean(log(c(1,0.33))),
769 mean(log(c(0.33,3)))))),length.out=20))
770 grid$z <- to.kcal(10^(grid$x*grid$y))
771 print(wireframe(z~x*y,grid,cuts=50,
774 scales=list(arrows=F,col=1,distance=0.75),
775 xlab=list(label=expression(paste(stdev,group("(",italic(log(cu[vesicle])),")"))),rot=30),
776 ylab=list(label=expression(paste(symbol("\341"),italic(log(cu[vesicle])),symbol("\361"))),rot=-35),
777 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
780 just=c("left","top"),
781 gp=gpar(fontsize=48),
782 y=unit(1,"npc")-unit(0.25,"lines"),
783 x=unit(0,"npc")+unit(0.25,"lines"))
786 \caption{Change in curvature forward ($cu_\mathrm{f}$) (A) or activation energy
787 in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the standard
788 deviation of the log of the curvature of the vesicle
789 ($\mathrm{stdev}\left(\log cu_\mathrm{vesicle}\right)$) and the mean
790 of the log of the curvature of the vesicle ($\left<\log
791 cu_\mathrm{vesicle}\right>$).}
792 \label{fig:cuf_graph}
795 \subsubsection{Length Forward}
797 As in the case of unsaturation, void formation is easier when vesicles
798 are made up of components of mismatched lengths. Thus, when the width
799 of the distribution of lengths is larger, the forward rate should be
800 greater as well, leading us to an equation of the form
801 $x^{\mathrm{stdev} l_\mathrm{vesicle}}$, where $\mathrm{stdev}
802 l_\mathrm{vesicle}$ is the standard deviation of the length of the
803 components of the vesicle, which has a maximum possible value of
804 $\Sexpr{format(digits=3,sd(c(rep(12,50),rep(24,50))))}$ and a minimum
805 of 0 in this set of simulations. Based on activation energy, a
806 reasonable base for $x$ is 2, leading to:
809 l_\mathrm{f} = 2^{\mathrm{stdev} l_\mathrm{vesicle}}
810 \label{eq:length_forward}
813 The most common $\mathrm{stdev} l_\mathrm{vesicle}$ is around $3.4$, which leads to
814 a $\Delta \Delta G^\ddagger$ of
815 $\Sexpr{format(digits=3,to.kcal(2^(3.4)))}
816 \frac{\mathrm{kcal}}{\mathrm{mol}}$.
818 % While it could be argued that increased length of the monomer could
819 % affect the rate of insertion into the membrane, it's not clear whether
820 % it would increase (by decreasing the number of available hydrogen
821 % bonds, for example) or decrease (by increasing the time taken to fully
822 % insert the acyl chain, for example) the rate of insertion or to what
823 % degree, so we do not take it into account in this formalism.
825 % \DLA{Incorporate McLean84 here}
826 % From McLean84LIB: Although it is difficult to measure cmc values for
827 % the sparingly soluble lipids used in this study, estimates for
828 % lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds,
829 % 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution;
830 % Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6
831 % X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X
832 % 10-8 M for DMPC was estimated from the linear relationship between ln
833 % cmc and the number of carbons in the PC acyl chain by using data for n
834 % = 7, 8, 10, and 16 [summarized in Tanford (1980)].
836 % From Nichols85: The magnitude of the dissociation rate constant
837 % decreases by a factor of approximately 3.2 per carbon increase in acyl
838 % chain length (see Table II here) {take into account for the formula;
841 % From Nichols85: The magnitude of the partition coefficient increases
842 % with acyl chain length [Keq(M-C6-NBD-PC) = (9.8±2.1) X 106 M-1 and Keq
843 % (P-C6-NBD-PC) = (9.4±1.3) x 107 M-1. {take into account for the
844 % formula; rz 8/17/2010}.
846 % From Nichols85: The association rate constant is independent of acyl
847 % chain length. {take into account for the formula; rz 8/17/2010}.
852 <<lf_graph,echo=FALSE,results='hide',fig.width=14,fig.height=5>>=
853 layout(matrix(1:2,nrow=1,ncol=2))
854 curve(2^x,from=0,to=sd(c(12,24)),
855 xlab=expression(paste(stdev,group("(",italic(l[vesicle]),")"))),
856 ylab=expression(italic(l[f])),
859 vps <- baseViewports()
860 pushViewport(vps$figure)
862 just=c("left","top"),
863 gp=gpar(fontsize=48),
864 y=unit(1,"npc")-unit(0.25,"lines"),
865 x=unit(0,"npc")+unit(0.25,"lines"))
867 curve(to.kcal(2^x),from=0,to=sd(c(12,24)),
868 xlab=expression(paste(stdev,group("(",italic(l[vesicle]),")"))),
869 ylab="Activation Energy (kcal/mol)",
871 vps <- baseViewports()
872 pushViewport(vps$figure)
874 just=c("left","top"),
875 gp=gpar(fontsize=48),
876 y=unit(1,"npc")-unit(0.25,"lines"),
877 x=unit(0,"npc")+unit(0.25,"lines"))
880 \caption{Change in length forward ($l_\mathrm{f}$) (A) or activation energy in
881 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the standard
882 deviation of the length of the vesicle
883 ($\mathrm{stdev}\left(l_\mathrm{vesicle}\right)$).}
887 \subsubsection{Complex Formation}
888 There is no contribution of complex formation to the forward reaction
889 rate in the current formalism.
893 \label{eq:complex_formation_forward}
896 \subsection{Backward adjustments ($k_{\mathrm{b}i\mathrm{adj}}$)}
898 Just as the forward rate constant adjustment $k_{\mathrm{f}i\mathrm{adj}}$
899 does, the backwards rate constant adjustment $k_{\mathrm{b}i\mathrm{adj}}$
900 takes into account unsaturation ($un_\mathrm{b}$), charge ($ch_\mathrm{b}$), curvature
901 ($cu_\mathrm{b}$), length ($l_\mathrm{b}$), and complex formation ($CF1_\mathrm{b}$), each of
902 which is modified depending on the specific component and the vesicle
903 from which the component is exiting:
907 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}
911 \subsubsection{Unsaturation Backward}
913 Unsaturation also influences the ability of a lipid molecule to leave
914 a membrane. If a molecule has an unsaturation level which is different
915 from the surrounding membrane, it will be more likely to leave the
916 membrane. The more different the unsaturation level is, the greater
917 the propensity for the lipid molecule to leave. However, a vesicle
918 with some unsaturation (eg. average unsaturation of 2) is more
919 favorable for lipids with more unsaturation (eg. 3) than lipids with
920 equivalently less unsaturation (eg. 1), so the difference in energy
921 between unsaturation is not linear. Therefore, an equation with the
922 shape $x^{\left| y^{-\left< un_\mathrm{vesicle}\right>
923 }-y^{-un_\mathrm{monomer}} \right| }$, where
924 $\left<un_\mathrm{vesicle}\right>$ is the average unsaturation of the
925 vesicle and $un_\mathrm{monomer}$ is the unsaturation of the monomer,
926 allows for increasing the efflux of molecules from membranes where
927 they strongly mismatch, while allowing vesicles with greater
928 unsaturation to tolerate greater unsaturation mismatch between
929 monomers and the vesicle. The bases x and y were chosen ad hoc to
930 produce reasonable Eyring energies of activation over the range of
931 unsaturations expected, leading to:
935 un_\mathrm{b} = 7^{1-\left(20\left(2^{-\left<un_\mathrm{vesicle} \right>} - 2^{-un_\mathrm{monomer}} \right)^2+1\right)^{-1}}
936 \label{eq:unsaturation_backward}
939 The most common $\left<un_\mathrm{vesicle}\right>$ is around $1.7$, which
940 leads to a range of $\Delta \Delta G^\ddagger$ from
941 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-0)^2+1))))}
942 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with 0 unsaturation
944 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-4)^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
945 for monomers with 4 unsaturations. See \cref{fig:unb_graph} for
946 the full range of possible values.
950 <<unb_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
951 trellis.device(new=F,color=TRUE)
952 trellis.par.set(list(axis.line =list(col="transparent")))
953 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
954 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
955 grid <- expand.grid(x=seq(0,4,length.out=20),
956 y=seq(0,4,length.out=20))
957 grid$z <- (7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1)))
958 print(wireframe(z~x*y,grid,cuts=50,
961 scales=list(arrows=F,col=1,distance=0.75),
962 xlab=list(label=expression(paste(symbol("\341"),italic(un[vesicle]),symbol("\361"))),rot=30),
963 ylab=list(label=expression(italic(un[monomer])),rot=-35),
964 zlab=list(label=expression(italic(un[b])),rot=93)),newpage=FALSE)
966 just=c("left","top"),
967 gp=gpar(fontsize=48),
968 y=unit(1,"npc")-unit(0.25,"lines"),
969 x=unit(0,"npc")+unit(0.25,"lines"))
971 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
972 grid <- expand.grid(x=seq(0,4,length.out=20),
973 y=seq(0,4,length.out=20))
974 grid$z <- to.kcal((7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1))))
975 print(wireframe(z~x*y,grid,cuts=50,
978 scales=list(arrows=F,col=1,distance=0.75),
979 xlab=list(label=expression(paste(symbol("\341"),italic(un[vesicle]),symbol("\361"))),rot=30),
980 ylab=list(label=expression(italic(un[monomer])),rot=-35),
981 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
984 just=c("left","top"),
985 gp=gpar(fontsize=48),
986 y=unit(1,"npc")-unit(0.25,"lines"),
987 x=unit(0,"npc")+unit(0.25,"lines"))
990 \caption{Change in unsaturation backward ($un_\mathrm{b}$) (A) or activation energy in
991 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the unsaturation of the
992 monomer leaving ($un_\mathrm{monomer}$) and the average unsaturation of
993 the vesicle ($\left<un_\mathrm{vesicle}\right>$).}
994 \label{fig:unb_graph}
997 \subsubsection{Charge Backwards}
998 As in the case of monomers entering a vesicle, opposites attract.
999 Monomers leaving a vesicle leave faster if their charge has the same
1000 sign as the average charge vesicle. An equation of the form
1001 $ch_\mathrm{b} = a^{\left<ch_v\right> ch_m}$ is then appropriate, and
1002 using a base of $a=20$ yields:
1005 ch_\mathrm{b} = 20^{\left<{ch}_v\right> {ch}_m}
1006 \label{eq:charge_backwards}
1009 The most common $\left<ch_v\right>$ is around $-0.164$, which leads to
1010 a range of $\Delta \Delta G^\ddagger$ from
1011 $\Sexpr{format(digits=3,to.kcal(20^(-.164*-1)))}
1012 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with charge $-1$ to
1013 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with charge $0$.
1014 See \cref{fig:chb_graph} for the full range of possible values
1019 <<chb_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1020 trellis.device(new=F,color=TRUE)
1021 trellis.par.set(list(axis.line =list(col="transparent")))
1022 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1023 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1024 x <- seq(-1,0,length.out=20)
1025 y <- seq(-1,0,length.out=20)
1026 grid <- expand.grid(x=x,y=y)
1027 grid$z <- as.vector(20^(outer(x,y)))
1028 print(wireframe(z~x*y,grid,cuts=50,
1031 scales=list(arrows=F,col=1,distance=0.75),
1032 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
1033 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
1034 zlab=list(label=expression(italic(ch[b])),rot=93)),newpage=FALSE)
1037 just=c("left","top"),
1038 gp=gpar(fontsize=48),
1039 y=unit(1,"npc")-unit(0.25,"lines"),
1040 x=unit(0,"npc")+unit(0.25,"lines"))
1042 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1043 x <- seq(-1,0,length.out=20)
1044 y <- seq(-1,0,length.out=20)
1045 grid <- expand.grid(x=x,y=y)
1046 grid$z <- to.kcal(as.vector(20^(outer(x,y))))
1047 print(wireframe(z~x*y,grid,cuts=50,
1050 scales=list(arrows=F,col=1,distance=0.75),
1051 xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30),
1052 ylab=list(label=expression(italic(ch[monomer])),rot=-35),
1053 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1056 just=c("left","top"),
1057 gp=gpar(fontsize=48),
1058 y=unit(1,"npc")-unit(0.25,"lines"),
1059 x=unit(0,"npc")+unit(0.25,"lines"))
1062 \caption{Change in charge backward ($ch_\mathrm{b}$) (A) or activation energy in
1063 $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the charge of the
1064 monomer leaving ($ch_\mathrm{monomer}$) and the average charge of
1065 the vesicle ($\left<ch_\mathrm{vesicle}\right>$).}
1066 \label{fig:chb_graph}
1070 \subsubsection{Curvature Backwards}
1072 The less a monomer's intrinsic curvature matches the average curvature
1073 of the vesicle in which it is in, the greater its rate of efflux. If
1074 the curvatures match exactly, $cu_\mathrm{f}$ needs to be one. To map negative and
1075 positive curvature to the same range, we also need take the logarithm.
1076 Positive ($cu > 1$) and negative ($0 < cu < 1$) curvature lipids cancel each other out,
1077 resulting in an average curvature of 1; the average of the log also
1078 has this property (average curvature of 0). Increasing mismatches in curvature increase the
1079 rate of efflux, but asymptotically. An equation which satisfies these
1080 criteria has the form $cu_\mathrm{f} = a^{1-\left(b\left( \left< \log
1081 cu_\mathrm{vesicle} \right> -\log
1082 cu_\mathrm{monomer}\right)^2+1\right)^{-1}}$. An alternative
1083 form would use the absolute value of the difference, however, this
1084 yields a cusp and sharp increase about the curvature equilibrium. We
1085 have chosen bases of $a=7$ and $b=20$ ad hoc, based on activation
1086 energy considerations.
1089 cu_\mathrm{b} = 7^{1-\left(20\left(\left< \log cu_\mathrm{vesicle} \right> -\log cu_\mathrm{monomer} \right)^2+1\right)^{-1}}
1090 \label{eq:curvature_backwards}
1093 The most common $\left<\log cu_\mathrm{vesicle}\right>$ is around
1094 $-0.013$, which leads to a range of $\Delta \Delta G^\ddagger$ from
1095 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(0.8))^2+1))))}
1096 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature 0.8 to
1097 $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(1.3))^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
1098 for monomers with curvature 1.3 to
1099 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature near
1100 1. The full range of values possible for $cu_\mathrm{b}$ are shown in
1101 \cref{fig:cub_graph}.
1103 % \RZ{What about the opposite curvatures that actually do fit to each
1107 % From Kumar91: It has also been shown experimentally that molecules
1108 % having complementary molecular shapes can form bilayer structures.
1109 % Inverted cone (wedge)- shaped lysoPC and cone-shaped Chol or
1110 % unsaturated PE can interact stoichiometrically to form bilayer
1111 % structures (9-12). {quote to justify our formula for complementary
1116 <<cub_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1117 trellis.device(new=F,color=TRUE)
1118 trellis.par.set(list(axis.line =list(col="transparent")))
1119 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1120 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1121 grid <- expand.grid(x=log(seq(0.8,1.33,length.out=20)),
1122 y=seq(0.8,1.33,length.out=20))
1123 grid$z <- 7^(1-1/(20*(grid$x-log(grid$y))^2+1))
1124 print(wireframe(z~x*y,grid,cuts=50,
1127 scales=list(arrows=F,col=1,distance=0.75),
1128 xlab=list(label=expression(paste(symbol("\341"),log(italic(cu[vesicle])),symbol("\361"))),rot=30),
1129 ylab=list(label=expression(italic(cu[monomer])),rot=-35),
1130 zlab=list(label=expression(italic(cu[b])),rot=93)),newpage=FALSE)
1133 just=c("left","top"),
1134 gp=gpar(fontsize=48),
1135 y=unit(1,"npc")-unit(0.25,"lines"),
1136 x=unit(0,"npc")+unit(0.25,"lines"))
1138 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1139 grid <- expand.grid(x=log(seq(0.8,1.33,length.out=20)),
1140 y=seq(0.8,1.33,length.out=20))
1141 grid$z <- to.kcal(7^(1-1/(20*(grid$x-log(grid$y))^2+1)))
1142 print(wireframe(z~x*y,grid,cuts=50,
1145 scales=list(arrows=F,col=1,distance=0.75),
1146 xlab=list(label=expression(paste(symbol("\341"),log(italic(cu[vesicle])),symbol("\361"))),rot=30),
1147 ylab=list(label=expression(italic(cu[monomer])),rot=-35),
1148 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1151 just=c("left","top"),
1152 gp=gpar(fontsize=48),
1153 y=unit(1,"npc")-unit(0.25,"lines"),
1154 x=unit(0,"npc")+unit(0.25,"lines"))
1157 \caption{Change in curvature backward ($cu_\mathrm{b}$) (A) or activation
1158 energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the
1159 curvature of the monomer leaving ($ch_\mathrm{monomer}$) and the
1160 average of the log of the curvature of the vesicle
1161 ($\left<\log cu_\mathrm{vesicle}\right>$).}
1162 \label{fig:cub_graph}
1166 \subsubsection{Length Backwards}
1168 In a model membrane, the dissociation constant increases by a factor
1169 of approximately 3.2 per carbon decrease in acyl chain length~%
1170 \citep{Nichols1985:phospholipid_monomer_vesicle_thermodynamics}.
1171 Unfortunately, the experimental data known to us only measures the
1172 effect of chain lengths less than or equal to the bulk lipid, not the
1173 effect of lengths exceeding it, and is only known for one bulk lipid
1174 component (DOPC). We assume then, that the increase is in relationship
1175 to the average vesicle, and that lipids with larger acyl chain length
1176 will also show an increase in their dissociation constant.
1179 l_\mathrm{b} = 3.2^{\left|\left<l_\mathrm{vesicle}\right>-l_\mathrm{monomer}\right|}
1180 \label{eq:length_backward}
1183 The most common $\left<l_\mathrm{vesicle}\right>$ is around $17.75$,
1184 which leads to a range of $\Delta \Delta G^\ddagger$ from
1185 $\Sexpr{format(digits=3,to.kcal(3.2^abs(12-17.75)))}
1186 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with length 12
1187 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$
1188 for monomers with length near 18 to
1189 $\Sexpr{format(digits=3,to.kcal(3.2^abs(24-17.75)))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
1190 for monomers with length 24. The full range of possible values of
1191 $l_\mathrm{b}$ are shown in \cref{fig:lb_graph}
1193 % (for methods? From McLean84LIB: The activation free energies and free
1194 % energies of transfer from self-micelles to water increase by 2.2 and
1195 % 2.1 kJ mol-' per methylene group, respectively. \RZ{see if we can use it
1196 % to justify arranging our changed in activating energy around 1
1201 <<lb_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1202 trellis.device(new=F,color=TRUE)
1203 trellis.par.set(list(axis.line =list(col="transparent")))
1204 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1205 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1206 grid <- expand.grid(x=seq(12,24,length.out=20),
1207 y=seq(12,24,length.out=20))
1208 grid$z <- 3.2^(abs(grid$x-grid$y))/1e6
1209 print(wireframe(z~x*y,grid,cuts=50,
1212 scales=list(arrows=F,col=1,distance=0.75),
1213 xlab=list(label=expression(paste(symbol("\341"),italic(l[vesicle]),symbol("\361"))),rot=30),
1214 ylab=list(label=expression(italic(l[monomer])),rot=-35),
1215 zlab=list(label=expression(italic(l[b])%*%10^-6),rot=93)),newpage=FALSE)
1218 just=c("left","top"),
1219 gp=gpar(fontsize=48),
1220 y=unit(1,"npc")-unit(0.25,"lines"),
1221 x=unit(0,"npc")+unit(0.25,"lines"))
1223 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1224 grid <- expand.grid(x=seq(12,24,length.out=20),
1225 y=seq(12,24,length.out=20))
1226 grid$z <- to.kcal(3.2^(abs(grid$x-grid$y)))
1227 print(wireframe(z~x*y,grid,cuts=50,
1230 scales=list(arrows=F,col=1,distance=0.75),
1231 xlab=list(label=expression(paste(symbol("\341"),italic(l[vesicle]),symbol("\361"))),rot=30),
1232 ylab=list(label=expression(italic(l[monomer])),rot=-35),
1233 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1236 just=c("left","top"),
1237 gp=gpar(fontsize=48),
1238 y=unit(1,"npc")-unit(0.25,"lines"),
1239 x=unit(0,"npc")+unit(0.25,"lines"))
1242 \caption{Change in length backwards ($l_\mathrm{b}\times 10^{-6}$) (A) or activation energy
1243 in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the length of the
1244 monomer leaving ($l_\mathrm{monomer}$) and the average length of
1245 lipids in the vesicle ($\left<l_\mathrm{vesicle}\right>$).}
1246 \label{fig:lb_graph}
1250 \subsubsection{Complex Formation Backward}
1252 Complex formation ($CF1$) describes the interaction between \ac{CHOL} and
1253 \ac{PC} or \ac{SM}, where \ac{PC} or \ac{SM} protects the hydroxyl group of \ac{CHOL} from
1254 interactions with water %
1255 \citep{Huang1999:chol_solubility_model,
1256 Huang1999:cholesterol_solubility, McConnell2006, McConnell2003}%
1257 . \ac{PC} ($CF1=2$) can interact with two \ac{CHOL}, and \ac{SM} ($CF1=3$) with three
1258 \ac{CHOL} ($CF1=-1$). If the average of $CF1$ is positive (excess of \ac{SM} and
1259 \ac{PC} with regards to complex formation), components with negative $CF1$
1260 (\ac{CHOL}) will be retained. If average $CF1$ is negative, components with
1261 positive $CF1$ are retained. An equation which has this property is
1262 $CF1_\mathrm{b}=a^{\left<CF1_\mathrm{vesicle}\right>
1263 CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{vesicle}\right>
1264 CF1_\mathrm{monomer}\right|}$, where difference of the exponent is
1265 zero if the average $CF1$ and the $CF1$ of the component have the same
1266 sign, or double the product if the signs are different. Based on
1267 activation energy considerations, we took an ad hoc base for $a$ as
1272 CF1_\mathrm{b}=1.5^{\left<CF1_\mathrm{vesicle}\right> CF1_\mathrm{monomer}-\left|\left<CF1_\mathrm{vesicle}\right> CF1_\mathrm{monomer}\right|}
1273 \label{eq:complex_formation_backward}
1276 The most common $\left<CF1_\mathrm{vesicle}\right>$ is around $0.925$,
1277 which leads to a range of $\Delta \Delta G^\ddagger$ from
1278 $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*-1-abs(0.925*-1))))}
1279 \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
1281 $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*2-abs(0.925*2))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$
1282 for monomers with complex formation $2$ to
1283 $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex
1284 formation $0$. The full range of possible values for $CF1_\mathrm{b}$ are
1285 depicted in \cref{fig:cf1b_graph}.
1290 <<cf1b_graph,echo=FALSE,results='hide',fig.width=14,fig.height=7>>=
1291 trellis.device(new=F,color=TRUE)
1292 trellis.par.set(list(axis.line =list(col="transparent")))
1293 pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off"))
1294 pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off"))
1295 grid <- expand.grid(x=seq(-1,3,length.out=20),
1296 y=seq(-1,3,length.out=20))
1297 grid$z <- 1.5^(grid$x*grid$y-abs(grid$x*grid$y))
1298 print(wireframe(z~x*y,grid,cuts=50,
1301 scales=list(arrows=F,col=1,distance=0.75),
1302 xlab=list(label=expression(paste(symbol("\341"),italic(CF1[vesicle]),symbol("\361"))),rot=30),
1303 ylab=list(label=expression(italic(CF1[monomer])),rot=-35),
1304 zlab=list(label=expression(italic(CF1[b])),rot=93)),newpage=FALSE)
1307 just=c("left","top"),
1308 gp=gpar(fontsize=48),
1309 y=unit(1,"npc")-unit(0.25,"lines"),
1310 x=unit(0,"npc")+unit(0.25,"lines"))
1312 pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off"))
1313 grid <- expand.grid(x=seq(-1,3,length.out=20),
1314 y=seq(-1,3,length.out=20))
1315 grid$z <- to.kcal(1.5^(grid$x*grid$y-abs(grid$x*grid$y)))
1316 print(wireframe(z~x*y,grid,cuts=50,
1319 scales=list(arrows=F,col=1,distance=0.75),
1320 xlab=list(label=expression(paste(symbol("\341"),italic(CF1[vesicle]),symbol("\361"))),rot=30),
1321 ylab=list(label=expression(italic(CF1[monomer])),rot=-35),
1322 zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE)
1325 just=c("left","top"),
1326 gp=gpar(fontsize=48),
1327 y=unit(1,"npc")-unit(0.25,"lines"),
1328 x=unit(0,"npc")+unit(0.25,"lines"))
1331 \caption{Change in complex formation 1 backwards ($CF1_\mathrm{b}$) (A) or
1332 activation energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus
1333 the complex formation of the monomer leaving
1334 ($CF1_\mathrm{monomer}$) and the average complex formation of the
1335 vesicle ($\left<CF1_\mathrm{vesicle}\right>$).}
1336 \label{fig:cf1b_graph}
1339 \section{Simulation Methodology}
1341 \subsection{Overall Architecture}
1343 The simulations are currently run using a program written in perl with
1344 various modules to handle the subsidiary parts. It produces output for
1345 each generation, including the step immediately preceding and
1346 immediately following a vesicle split (and optionally, each step) that
1347 is written to a state file which contains the state of the vesicle,
1348 environment, kinetic parameters, program invocation options, time, and
1349 various other parameters necessary to recreate the state vector at a
1350 given time. This output file is then read by a separate program in
1351 perl to produce different output which is generated out-of-band;
1352 output which includes graphs and statistical analysis is performed
1353 using R~\citep{R:maincite} (and various grid graphics modules) called
1354 from the perl program.
1356 The separation of simulation and output generation allows refining
1357 output, and simulations performed on different versions of the
1358 underlying code to be compared using the same analysis methods and
1359 code. It also allows later simulations to be restarted from a specific
1360 generation, utilizing the same environment.
1362 Random variables of different distributions are calculated using
1363 Math::Random (0.71), which is seeded for each run using entropy from
1364 the linux kernel's urandom device.
1366 Code is available upon request.
1368 \subsection{Environment Creation}
1370 \subsubsection{Components}
1371 The environment contains concentrations of components. In the current
1372 set of simulations, there are \Sexpr{1+4*5*7} different possible
1373 components, consisting of \ac{PC}, PE, \ac{PS}, \ac{SM}, and \ac{CHOL}, with all lipids
1374 except for \ac{CHOL} having 5 possible unsaturations ranging from 0 to 4,
1375 and 7 lengths from $12,14,...,22$ ($7\cdot
1376 5\cdot4+1=\Sexpr{1+4*5*7}$). In cases where the environment has less
1377 than the maximum number of components, the components are selected in
1378 order without replacement from a randomly shuffled deck of components
1379 (with the exception of \ac{CHOL}) represented once until the desired number
1380 of unique components is obtained. \ac{CHOL} is overrepresented
1381 $\Sexpr{5*7}$ times to be at the level of other lipid types, ensuring
1382 that the probability of \ac{CHOL} being absent in the environment is the
1383 same as the probability of any one of the other lipid types (\ac{PS}, PE,
1384 etc.) being absent. This reduces the number of simulations with a
1385 small number of components which are devoid of \ac{CHOL}.
1387 \subsubsection{Concentrations}
1388 Once the components of the environment have been selected, their
1389 concentrations are determined. In experiments where the environmental
1390 concentration is the same across all lipid components, the
1391 concentration is set to $10^{-10}$~M. In other cases, the
1392 environmental concentration is set to a random number from a gamma
1393 distribution with shape parameter 2 and an average of
1394 $10^{-10}$~M. The base concentration ($10^{-10}$~M)
1395 can also be altered at the initialization of the experiment to
1396 specific values for specific lipid types or components.
1398 The environment is a volume which is the maximum number of vesicles
1399 from a single simulation (4096) times the size of the vesicle
1400 (usually 10000) divided by Avagadro's number divided by the total
1401 environmental lipid concentration, or usually
1402 \Sexpr{4096*10000/6.022E23/141E-10}~L.
1404 \subsection{Initial Vesicle Creation}
1406 Initial vesicles are seeded by selecting lipid molecules from the
1407 environment until the vesicle reaches a specific starting size. The
1408 vesicle starting size has gamma distribution with shape parameter 2
1409 and a mean of the per-simulation specified starting size, with a
1410 minimum of 5 lipid molecules, or can be specified to have a precise
1411 number of molecules. Lipid molecules are then selected to be added to
1412 the lipid membrane according to four different methods. In the
1413 constant method, molecules are added in direct proportion to their
1414 concentration in the environment. The uniform method adds molecules in
1415 proportion to their concentration in the environment scaled by a
1416 uniform random value, whereas the random method adds molecules in
1417 proportion to a uniform random value. The final method is a binomial
1418 method, which adjust the probability of adding a molecule of a
1419 specific component by the concentration of that component, and then
1420 adds the molecules one by one to the membrane. This final method is
1421 also used in the first three methods to add any missing molecules to
1422 the starting vesicle which were unallocated due to the requirement to
1423 add integer numbers of molecules. (For example, if a vesicle was to
1424 have 10 molecules split evenly between three lipid types, the 10th
1425 molecule would be assigned by randomly choosing between the three
1426 lipid types, weighted by their concentration in the environment.)
1428 \subsection{Simulation Step}
1430 Once the environment has been created and the initial vesicle has been
1431 formed, molecules join and leave the vesicle based on the kinetic
1432 parameters and state equation discussed above until the vesicle splits
1433 forming two progeny vesicles. The program then follows both vesicles
1434 and their progeny until the simulation reaches either the maximum
1435 number of vesicles (usually 4096), or the maximum simulation time
1438 \subsubsection{Calculation of Vesicle Properties}
1440 \label{sec:ves_prop}
1441 $S_\mathrm{vesicle}$ is the surface area of the vesicle, and is the sum of
1442 the surface area of all of the individual lipid components; each lipid
1443 type has a different surface area; we assume that the lipid packing is
1444 optimal, and there is no empty space.
1446 For the other vesicle properties (length, unsaturation, charge, and
1447 curvature), we calculate the mean, the standard deviation, the mean of
1448 the log, the mean of the absolute value of the log, the standard
1449 deviation of the log, and the standard deviation of the absolute value
1450 of the log. All cases, when we take the log, we take the log of the
1451 absolute value, and map $\log(0)$ to $0$. For the purpose of plotting
1452 vesicle properties, we only plot the mean of the property.
1454 \subsubsection{Joining and Leaving of Lipid Molecules}
1456 Determining the number of molecules to add to the lipid membrane
1457 ($n_i$) requires knowing $k_{\mathrm{f}i_\mathrm{adj}}$, the surface area of the
1458 vesicle $S_\mathrm{vesicle}$ (see \cref{sec:ves_prop}), the time interval
1459 $dt$ during which lipids are added, the base $k_{\mathrm{f}i}$, and the
1460 concentration of the monomer in the environment
1461 $\left[C_{i\mathrm{monomer}}\right]$ (see \cref{eq:state}).
1462 $k_{\mathrm{f}i\mathrm{adj}}$ is calculated (see \cref{eq:kf_adj}) based on the
1463 vesicle properties and their hypothesized effect on the rate (in as
1464 many cases as possible, experimentally based)
1465 (see \cref{sec:kinetic_adjustments} for details). $dt$ can be varied
1466 (see \cref{sec:step_duration}), but for a given step is constant. This
1467 leads to the following:
1469 $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$
1471 In the cases where $n_i > 1$, the integer number of molecules is
1472 added. Fractional $n_i$ or the fractional remainder after the addition
1473 of the integer molecules are the probability of adding a specific
1474 molecule, and are compared to a uniformly distributed random value
1475 between 0 and 1. If the random value is less than or equal to the
1476 fractional part of $n_i$, an additional molecule is added.
1478 Molecules leaving the vesicle are handled in a similar manner, with
1480 $n_i = k_{\mathrm{b}i}k_{\mathrm{b}i_\mathrm{adj}}C_{i_\mathrm{vesicle}}\mathrm{N_A}dt$.
1482 Before addition or removal of molecules, the properties of the vesicle
1483 are calculated; this is done so that molecules can be considered to be
1484 added and removed at the same instant, even though additions are
1485 handled first programmatically. This also avoids cases where a removal
1486 would have resulted in a negative number of molecules for a particular
1489 \subsubsection{Step duration}
1490 \label{sec:step_duration}
1492 $dt$ is the time taken for each step of joining, leaving, and checking
1493 split. In the current implementation, it starts with a value of
1494 $10^{-6}$~$\mathrm{s}$ but this is modified in between each step if the
1495 number of molecules joining or leaving is too large or small. If more
1496 than half of the starting vesicle size molecules join or leave in a
1497 single step, $dt$ is reduced by half. If less than the starting
1498 vesicle size molecules join or leave in 100 steps, $dt$ is doubled.
1499 This is necessary to curtail run times and to automatically adjust to
1500 different experimental runs.
1502 \subsubsection{Vesicle splitting}
1504 If a vesicle has grown to a size which is more than double the
1505 starting vesicle size, the vesicle splits. Molecules are assigned to
1506 the progeny vesicles at random, with each progeny vesicle having an
1507 equal chance of getting a single molecule. The number of molecules to
1508 assign to the first vesicle has binomial distribution with a
1509 probability of an event in each trial of 0.5 and a number of trials
1510 equal to the number of molecules.
1514 The environment, initial vesicle, and the state of the vesicle
1515 immediately before and immediately after splitting are stored
1516 to produce later output.
1518 % \section{Analyzing output}
1520 % The analysis of output is handled by a separate perl program which
1521 % shares many common modules with the simulation program. Current output
1522 % includes simulation progress, summary tables, summary statistics, and
1525 % \subsection{PCA plots}
1527 % Two major groups of axes that contribute to vesicle variation between
1528 % subsequent generations are the components and properties of vesicles.
1529 % Each component in a vesicle is an axis on its own; it can be measured
1530 % either as an absolute number of molecules in each component, or the
1531 % fraction of molecules of that component over the total number of
1532 % molecules; the second approach is often more convenient, as it allows
1533 % vesicles with different number of molecules to be directly compared
1534 % (though it hides any effect of vesicle size).
1536 % In order to visualize the variation between and transition of
1537 % subsequent generations of vesicles from their initial state in the
1538 % simulation, to their final state at the termination of the simulation,
1539 % we plot the projection of the generations onto the two principle PCA
1540 % axes (and alternatively, any pairing of the axes).
1542 % \subsection{Carpet plots}
1544 % Carpet plots show the distance/similarity between the composition or
1545 % properties of all generations in a single run. The difference between
1546 % each group of vesicle can be calculated using Euclidean distance
1547 % (properties and compositions) or H similarity (composition only). We
1548 % must use Euclidean distance for properties because the H distance
1549 % metric is invalid when comparing negative vector elements to positive
1552 % In addition to showing distance or similarity, carpet plots also
1553 % depict propersomes and composomes as square boxes on the diagonals and
1554 % propertypes and compotypes as rectangles off the diagonals, each
1555 % propertype or compotype with a distinct color.
1557 % \subsection{Propersomes, propertypes, composomes and compotypes}
1559 % A generation is considered to be a propersome if it is less than $0.6$
1560 % units (by Euclidean distance of normalized properties) away from the
1561 % generation immediately following and preceding. Likewise, a generation
1562 % is in a composome if its H similarity is more than $0.9$ (by the
1563 % normalized H metric) from the preceding and following generations.
1564 % Propersomes and composomes are then assigned to propertypes and
1565 % compotypes using Paritioning Around Medioids
1566 % (PAM). All values of $k$ between 2 and 15
1567 % (or the number of propersomes and composomes, if that is less than 15)
1568 % are tried, and the clustering with the smallest
1569 % silhouette~\citep{Rousseeuw1987:silhouettes} is chosen as the ideal
1570 % clustering~\citep{Shenhav2005:pgard}.
1573 %\bibliographystyle{unsrtnat}
1574 %\bibliography{references.bib}