\documentclass[english,12pt]{article} \usepackage{fancyhdr} %\usepackage[pdftex]{graphicx} \usepackage{graphicx} \usepackage[bf]{caption2} \usepackage{rotating} \usepackage{multirow} \usepackage{textcomp} \usepackage{mathrsfs} \usepackage{amssymb} \usepackage{setspace} \usepackage{txfonts} \usepackage[light,all]{draftcopy} \usepackage{fancyref} \usepackage{acronym} \usepackage{array} \usepackage{dcolumn} \usepackage{booktabs} \usepackage[noblocks]{authblk} \usepackage[hyperfigures,bookmarks,colorlinks,citecolor=black,filecolor=black,linkcolor=black,urlcolor=black]{hyperref} \usepackage[capitalise]{cleveref} \usepackage[sectionbib,sort&compress,numbers]{natbib} \usepackage[nomargin,inline,draft]{fixme} \usepackage[x11names,svgnames]{xcolor} \usepackage{texshade} \newenvironment{narrow}[2]{% \begin{list}{}{% \setlength{\topsep}{0pt}% \setlength{\leftmargin}{#1}% \setlength{\rightmargin}{#2}% \setlength{\listparindent}{\parindent}% \setlength{\itemindent}{\parindent}% \setlength{\parsep}{\parskip}}% \item[]}{\end{list}} \newenvironment{paperquote}{% \begin{quote}% \it }% {\end{quote}} \renewcommand{\textfraction}{0.15} \renewcommand{\topfraction}{0.85} \renewcommand{\bottomfraction}{0.65} \renewcommand{\floatpagefraction}{0.60} %\renewcommand{\baselinestretch}{1.8} \newenvironment{enumerate*}% {\begin{enumerate}% \setlength{\itemsep}{0pt}% \setlength{\parskip}{0pt}}% {\end{enumerate}} \newenvironment{itemize*}% {\begin{itemize}% \setlength{\itemsep}{0pt}% \setlength{\parskip}{0pt}}% {\end{itemize}} \newcommand{\DLA}[1]{\textcolor{red}{\fxnote{DLA: #1}}} \newcommand{\RZ}[1]{\textcolor{red}{\fxnote{RZ: #1}}} \newcommand{\OM}[1]{\textcolor{red}{\fxnote{OM: #1}}} \oddsidemargin 0.0in \textwidth 6.5in \raggedbottom \clubpenalty = 10000 \widowpenalty = 10000 \pagestyle{fancy} \author{} \title{Supplemental Material} \date{} % rubber: module bibtex % rubber: module natbib \onehalfspacing \begin{document} \maketitle <>= opts_chunk$set(dev="CairoPDF",out.width="\\columnwidth",out.height="0.7\\textheight",out.extra="keepaspectratio") opts_chunk$set(cache=TRUE, autodep=TRUE) options(scipen = -2, digits = 1) library("lattice") library("grid") library("Hmisc") library("gridBase") to.latex <- function(x){ gsub("\\\\","\\\\\\\\",latexSN(x)) } # R in cal / mol K to.kcal <- function(k,temp=300) { gasconst <- 1.985 return(-gasconst*temp*log(k)/1000) } @ % \setcounter{figure}{0} \setcounter{table}{0} \makeatletter \renewcommand{\thefigure}{S\@arabic\c@figure} \renewcommand{\thetable}{S\@arabic\c@table} \makeatother % \section{Competition Implementation} % \subsection{Implementation changes} % % \begin{itemize} % \item settable maximum number of vesicles to track (default $10^4$) % \item start with 1~L ($10^{-3}$~m$^3$) cube % \item if at any point the number of vesicles exceeds the maximum % number, chop the volume and environment molecule number into tenths, % randomly select one tenth of the vesicles, and continue tracking. % \item generations will be counted per vesicle, and each progeny % vesicle will have a generation number one greater than the parental % vesicle. % \item 100 generations can result in as many as $2^{100}$ % ($\Sexpr{2^100}$) vesicles or as few as % 101 vesicles. % \item Environment will use a specific number of each component instead % of a constant concentration; as the number may be larger than % \texttt{long long} ($2^{64}$), we use libgmp to handle an arbitrary % precision number of components % \end{itemize} % % \subsection{Infrastructure changes} % \begin{itemize} % \item Rewrite core bits in C % \item Use libgmp for handling large ints % \item Use openmpi to split the calculations out over multiple % machines/processors and allow deploying to larger % clusters/supercomputers % \end{itemize} \section{State Equation} % double check this with the bits in the paper For a system with monomers $(_\mathrm{monomer})$ and a vesicle $(_\mathrm{vesicle})$, the change in composition of the $i$th component of a lipid vesicle per change in time ($d \left[C_{i_\mathrm{vesicle}}\right]/dt$) can be described by a modification of the basic mass action law: \begin{equation} \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] - k_{\mathrm{b}i}k_{\mathrm{b}i\mathrm{adj}}\left[C_{i_\mathrm{vesicle}}\right] \label{eq:state} \end{equation} The base forward kinetic parameter for the $i$th component is $k_{\mathrm{f}i}$ and is dependent on the particular lipid type (PC, PS, SM, etc.). The forward adjustment parameter, $k_{\mathrm{f}i\mathrm{adj}}$, is based on the properties of the vesicle and the specific component (type, length, unsaturation, etc.) (see \cref{eq:kf_adj,sec:kinetic_adjustments}). $\left[C_{i_\mathrm{monomer}}\right]$ is the molar concentration of monomer of the $i$th component. $\left[S_\mathrm{vesicle}\right]$ is the surface area of the vesicle per volume. The base backwards kinetic parameter for the $i$th component is $k_{\mathrm{b}i}$ and its adjustment parameter $k_{\mathrm{b}i\mathrm{adj}}$ (see \cref{eq:kb_adj,sec:kinetic_adjustments}). $\left[C_{i_\mathrm{vesicle}}\right]$ is the molar concentration of the $i$th component in the vesicle. \subsection{Per-Lipid Kinetic Parameters} <>= kf.prime <- c(3.7e6,3.7e6,5.1e7,3.7e6,2.3e6) kf <- (as.numeric(kf.prime)*10^-3)/(63e-20*6.022e23) @ Each of the 5 lipid types has different kinetic parameters; where available, these were taken from literature. % \begin{table} % \centering % \begin{tabular}{c c c c c c c c} % \toprule % 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 \\ % \midrule % 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 \\ % 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 \\ % 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 \\ % 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 \\ % 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 \\ % \bottomrule % \end{tabular} % \caption{Kinetic parameters and molecular properties of lipid types} % \label{tab:kinetic_parameters_lipid_types} % \end{table} %%% \DLA{I think we may just reduce these three sections; area, $k_\mathrm{f}$ %%% and $k_\mathrm{b}$ to \cref{tab:kinetic_parameters_lipid_types} with %%% references.} %%% %%% \RZ{Yes, but we also have to have then as comments the numbers that %%% are not supported by refs} \subsubsection{$k_\mathrm{f}$ for lipid types} $k_{\mathrm{f}_\mathrm{PC}}$ was measured by \citet{Nichols1985:phospholipid_monomer_vesicle_thermodynamics} and was found to be $3.7\times 10^6$~$\frac{1}{\mathrm{M}\mathrm{s}}$ by the partitioning of P-C$_6$-NBD-PC between DOPC vesicles and water. As similar references are not available for SM or PS, we assume that they have the same $k_\mathrm{f}$. For CHOL, no direct measurement of $k_\mathrm{f}$ is available, however, \citet{Estronca2007:dhe_kinetics} measured the transfer of DHE from BSA to LUV, and found a $k_\mathrm{f}$ of $5.1\times 10^7$~% $\frac{1}{\mathrm{M} \mathrm{s}}$. We assume that this value is close to that of CHOL, and use it for $k_{\mathrm{f}_\mathrm{CHOL}}$. In the case of PE, \citet{Abreu2004:kinetics_ld_lo} measured the association of NBD-DMPE with POPC LUV found a value for $k_\mathrm{f}$ of $2.3 \times 10^{6}$~% $\frac{1}{\mathrm{M} \mathrm{s}}$. These three authors used a slightly different kinetic formalism ($\frac{d\left[A_v\right]}{dt} = k'_\mathrm{f}[A_m][L_v] - k_\mathrm{b}[A_v]$), so we correct their values of $k_\mathrm{f}$ by multiplying by the surface area of a mole of lipids. \subsubsection{$k_\mathrm{b}$ for lipid types} \citet{Wimley1990:dmpc_exchange} measured the half time for the exchange of [$^3$H]DMPC between LUV at 30\textdegree C, and found it to be 9.6 hr. As this is a first order reaction, and the primary limiting step in exchange is lipid desorption, $k_\mathrm{b}$ for DMPC is $k_{\mathrm{b}_\mathrm{PC}}=\frac{\log 2}{9.6 \times 60 \times 60} \approx \Sexpr{log(2)/(9.6*60*60)} \mathrm{s}^{-1}$. We assume that $k_\mathrm{b}$ for SM is the same as for PC. To estimate the $k_\mathrm{b}$ of PE and PS, we used the data from \citet{Nichols1982:ret_amphiphile_transfer} who measured the rate of exchange of the fluorescent label C$_6$-NBD attached to different lipid species. Although the values of kb are different for the labeled and unlabeled lipids, we assume that the ratios of the kinetics constants for the lipid types are the same. Furthermore we assume that PG behaves similarly to PS. Thus, we can determine the $k_\mathrm{b}$ of PE and PS from the already known $k_\mathrm{b}$ of PC. For C$_6$-NBD labeled PC, \citet{Nichols1982:ret_amphiphile_transfer} obtained a $k_\mathrm{b}$ of $0.89$~$\mathrm{min}^{-1}$, PE of $0.45$~$\mathrm{min}^{-1}$ and PG of $0.55$~$\mathrm{min}^{-1}$. Assuming the ratios of the rate of exchange are the same for unlabeled lipids and labeled lipids, we can determine the $k_\mathrm{b}$ of PE and PS from the already known $k_\mathrm{b}$ of PC~\citep{Wimley1990:dmpc_exchange}. Calculating the ratio, this leads us to $k_{\mathrm{b}_\mathrm{PE}} = \frac{k_{\mathrm{b}_\mathrm{PC}}\times\mathrm{PE}}{\mathrm{PC}} \approx \frac{2\times 10^{-5}\,\mathrm{s}^{-1} \times 0.45\,\mathrm{min}^{-1}}{0.89\,\mathrm{min}^{-1}} \approx \Sexpr{log(2)/(9.6*60*60)*0.45/0.89}$~$\mathrm{s}^{-1}$ and likewise, $k_{\mathrm{b}_\mathrm{PS}}\approx \Sexpr{log(2)/(9.6*60*60)*0.55/0.89}$~$\mathrm{s}^{-1}$. The $k_\mathrm{b}$ of SM was determined using the work of \citet{Bai1997:lipid_movementbodipy}, who measured spontaneous transfer of C$_5$-DMB-SM and C$_5$-DMB-PC from donor and acceptor vesicles, finding $3.4\times10^{-2}$~$\mathrm{s}^{-1}$ and $2.2\times10^{-3}$~$\mathrm{s}^{-1}$ respectively; using the ratio of $k_\mathrm{b}$ of C$_5$-DMB-SM to the $k_\mathrm{b}$ of C$_5$-DMB-PC times the $k_\mathrm{b}$ of PC ($\frac{3.4 \times 10^{-2}\mathrm{s}^{-1}}{2.2 \times 10^{-3}\mathrm{s}^{-1}}\approx \Sexpr{log(2)/(9.6*60*60)}\mathrm{s}^{-1}$), we obtain $k_{\mathrm{b}_\mathrm{SM}} \approx \Sexpr{log(2)/(9.6*60*60)*3.4e-2/2.2e-3}$. In the case of CHOL, \citet{Jones1990:lipid_transfer} measured the $t_{1/2}$ of [$^3$H] transfer from POPC vesicles and found it to be 41 minutes, leading to a $k_{\mathrm{b}_\mathrm{CHOL}} = \frac{\log 2}{41\times 60} \approx \Sexpr{log(2)/(41*60)}$~$\mathrm{s}^{-1}$. \subsubsection{Headgroup Surface Area for lipid types} % Smondyrev, A. M., and M. L. Berkowitz. 1999b. United atom force % field for phospholipid membranes: constant pressure molecular % dynamics simulation of dipalmitoylphosphatidicholine/water system. % J. Comput. Chem. 50:531–545. Different lipids have different headgroup surface areas, which contributes to $\left[S_\mathrm{vesicle}\right]$. \citet{Smaby1997:pc_area_with_chol} measured the surface area of POPC with a Langmuir film balance, and found it to be 63~Å$^2$ at $30$~$\frac{\mathrm{mN}}{\mathrm{m}}$. Molecular dynamic simulations found an area of 54 Å$^2$ for DPPS\citep{Cascales1996:mds_dpps_area,Pandit2002:mds_dpps}, which is in agreement with the experimental value of 56~Å$^2$ found using a Langmuir balance by \citet{Demel1987:ps_area}. \citet{Shaikh2002:pe_phase_sm_area} measured the area of SM using a Langmuir film balance, and found it to be 61~Å$^2$. Using $^2$H NMR, \citet{Thurmond1991:area_of_pc_pe_2hnmr} found the area of DPPE-d$_{62}$ to be 55.4 Å$^2$. \citet{Robinson1995:mds_chol_area} found an area for CHOL of 38~Å$^2$ using molecular dynamic simulations. % robinson's chol area is kind of crappy; they did it using MDS, but % vaguely; they don't indicate what it was at the end, and they only % had 2 molecules of CHOL. % PC: 63 A (Smaby97rd) (Pandit (?) % PS: 54 A (Pandit02LIB) (Cascales 1996; J. Chem. Phys 104:2713, Mukhopadhyay2004) % CHOL: 38 A (Robinson 1995) (Lundberg 1982) % SM: 51 A (Shaikh2002, but 61(?)) % PE: 55 A (Thurmond, 1991) (Pandit2002) % % Compare to results by Petrache2004 with DOPC of 72.5, DOPS of 65.3. % % The cross‑sectional area of cholesterol is ~37 ֵ2 (Lundberg, B. 1982. A % surface film study of the lateral packing of phosphatidylcholine and % cholesterol. Chem. Phys. Lipids. 31:23-32). % % From Pandit02LIB: the area per headgroup for PE is ~10‑20\% (Thurmond % et al., 1991) smaller than the area per PC molecule. For DLPE, it is % 50.6 ֵ2 (McIntosh and Simon, 1986; Damodaran et al., 1992). The % estimated area per molecule for DPPE is ~55.4 ֵ2 (Thurmond et al., % 1991). % % From Pandit02LIB: (from their simulations) average area per headgroup % for DPPS molecules is 53.75 ± 0.10 ֵ2. The values inferred from the % experiments are in the region between 45 and 55 ֵ2 (Cevc et al., 1981. % Titration of the phase transition of phosphatidylserine bilayer % membranes. Effect of pH, surface electrostatics, ion binding, and % head‑group hydration. Biochemistry. 20:4955‑4965; | Demel et al., % 1987. Monolayer characteristics and thermal behavior of natural and % synthetic phosphatidylserines. Biochemistry. 26:8659‑8665). Note that % the average area per headgroup for DPPC bilayer obtained from % simulations and measurements is ~62 ֵ A2. One would expect that the % area per headgroup in case of DPPS molecules will be larger than the % area per DPPC molecule, because the headgroups of DPPS are charged and % therefore repel each other. Contrary to this expectation, the area per % headgroup in DPPS is ~13\% smaller than that of DPPC. López Cascales % et al. (1996. Molecular dynamics simulation of a charged biological % membrane. J. Chem. Phys. 104:2713‑2720.) proposed that this reduction % in the area per headgroup is due to a strong intermolecular % coordination between DPPS molecules. MD gave area per POPS of 55 A2 at % 300K (Mukhopadhyay04LIB). % % However, by 2H NMR and X-ray, Petrache04LIB determined the area of % DOPS at 30 C to be 65.3 A2, considerably smaller than DOPC, determined % by the same group in another paper to be 72.5 ֵ A2 % (Tristram‑Nagle00LIB). % % From Shaikh02: At 35 mN/m of surface pressure, the areas/molecule for % POPE, SM, and CHOL were found to be 65.2, 61.1, and 38.8 Å2. % % In Lonchin99: POPC area 72 A2 ((a) Cornell, B. A.; Middlehurst, J.; % Separovic, F. Biochim. Biophys. Acta1980, 598, 405−410), oleic acid % area 32 A2 ((a) Small, D. M. In Handbook of Lipid Research; Hanahan, % D. J., Ed.; Plenum Press:  New York, 1986; Vol. 4, Chapter 9. % % (b) Jain, M. K. Introduction to Biological Membranes, 2nd ed.; John % Wiley \& Sons:  New York, 1988; Chapter 3). Cistola, D. P.; Atkinson, % D.; Hamilton, J. A.; Small, D. M. Biochemistry1986, 25, 2804−2812. % % From Marsh96: The values obtained for the headgroup area per molecule, % fitted as free parameters in Fig. 6, are APE = 58 A2 and APC = 78 A2, % for DOPE and DOPC, respectively. These values are within the range % obtained for the corresponding areas per molecule of % phosphatidylethanolamines and phosphatidylcholines in the fluid % lamellar phase (see, e.g., Marsh, 1990), and correspond to values for % the lipid packing parameter of V/Al = 1.08 and 1.34 for DOPC and DOPE, % respectively. The former value is consistent with DOPC alone forming % hydrated lamellar structures rather than nonlamellar structures. The % latter value lies at the lower end of the packing parameters obtained % for phosphatidylethanolamines. % % From Kumar91: The head group area chosen was 71.7 A2 for PC and % lysoPC, 42 A2 for PE, and 19 A2 for Chol (5, 14). The value selected % for PC is in agreement with the 71 A2 value determined for PCs with % saturated acyl chains (14). The hydrocarbon volume and the length % chosen for Chol are 400 A3 and 17.25 A, respectively (5). % % % % From Sampaio05: Besides this work and our own earlier report on the % association of NBD-DMPE with lipid bilayers (Abreu et al., 2004), we % are aware of only one other report in the literature (Nichols, 1985) % in which all the kinetic constants of lipid-derived amphiphile % association with lipid bilayer membranes were experimentally measured. % {indeed; everything is k- !!!; rz} % % From McLean84LIB: Although it is difficult to measure cmc values for % the sparingly soluble lipids used in this study, estimates for % lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds, % 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution; % Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6 % X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X % 10-8 M for DMPC was estimated from the linear relationship between ln % cmc and the number of carbons in the PC acyl chain by using data for n % = 7, 8, 10, and 16 [summarized in Tanford (1980)]. % % From Toyota08: Recently, several research groups have reported % self-reproducing systems of giant vesicles that undergo a series of % sequential transformations: autonomous growth, self-division, and % chemical reactions to produce membrane constituents within the giant % vesicles.44-47 % % Vesicle sizes of 25 nm for SUV and 150 nm for LUV were mentioned by % Thomas02. % % From Lund-Katz88: Charged and neutral small unilamellar vesicles % composed of either saturated PC, unsaturated PC, or SM had similar % size distributions with diameters of 23 \& 2 nm. % % From Sampaio05LIB: The exchange of lipids and lipid derivatives % between lipid bilayer vesicles has been studied for at least the last % 30 years. Most of this work has examined the exchange of amphiphilic % molecules between a donor and an acceptor population. The measured % efflux rates were shown in almost all cases, not surprisingly, to be % first order processes. In all of this work, the insertion rate has % been assumed to be much faster than the efflux rate. Having measured % both the insertion and desorption rate constants for amphiphile % association with membranes, our results show that this assumption is % valid. In several cases reported in the literature, the insertion rate % constant was assumed, although never demonstrated, to be a % diffusion-controlled process. % % (for methods? From McLean84LIB: The activation free energies and free % energies of transfer from self-micelles to water increase by 2.2 and % 2.1 kJ mol-' per methylene group, respectively. {see if we can use it % to justify arranging our changed in activating energy around 1 % kcal/mol; rz}). % % Jones90 give diameter of LUV as 100 nm, and of SUV as 20 nm; that % would give the number of molecules per outer leaflet of a vesicle as % 1500. % % Form Simard08: Parallel studies with SUV and LUV revealed that % although membrane curvature does have a small effect on the absolute % rates of FA transfer between vesicles, the ΔG of membrane desorption % is unchanged, suggesting that the physical chemical properties which % govern FA desorption are dependent on the dissociating molecule rather % than on membrane curvature. However, disagreements on this fundamental % issue continue (57, 61, 63, 64) % % (methods regarding the curvature effect: Kleinfeld93 showed that the % transfer parameters of long-chain FFA between the lipid vesicles % depend on vesicle curvature and composition. Transfer of stearic acid % is much slower from LUV as compared to SUV). % % From McLean84: In a well-defined experimental system consisting of % unilamellar lipid vesicles, in the absence of protein, the % rate-limiting step for the overall exchange process is desorption % (McLean \& Phillips, 1981). {thus I can take exchange data for the % estimation of k- rz; 8/11/08}. \subsubsection{Complex Formation 1} \citet{Thomas1988:chol_transfer} found that SM significantly decreases the rate of cholesterol transfer from PC, while PE and PS have no effect at physiologically significant levels. We attribute this to the formation of complexes between SM and CHOL, which retards the rate of efflux of CHOL from the membrane. Similar complexes exist between PC and CHOL, where it was shown that two CHOL molecules complex with one PC~\citep{Huang1999:chol_solubility_model, Huang1999:cholesterol_solubility,McConnell2006,McConnell2003}. It has also been shown that SM binds more CHOL molecules than PC~\citep{Katz1988:pl_packing_chol}. We assume that three CHOL complex with one SM, leading to $\mathrm{CF}1$ values of $3$, $2$, and $-1$ for SM, PC, and CHOL, respectively. \subsubsection{Curvature} We used the formulation for curvature of \citet{Israelachvili1975:amphiphile_self_assembly}, or $S=\frac{v}{l_c a}$, where $S$ is the curvature, which ranges from $0$ to $\infty$, $l_c$ is the critical length of the acyl chain, $v$ is the volume of the lipid, and $a$ is the area of the head group. \citet{Kumar1991:lipid_packing} found a value for $S$ of 1.33 for PE, 1.21 for CHOL, and 0.8 for diC$_{16}$ PC. We assume that PS has neutral curvature (value of 1), and that SM has the same curvature as PC (0.8). It should also be noted that in reality the curvature parameter changes with length, but at longer chain lengths, is effectively constant; in the current model, curvature is held constant for each species. % % 1. Israelachvili, J. N., Mitchell, D. J. & Ninham, B. W. (1976) J. % Chem. Soc. Faraday Trans. 2 72, 1525-1568. % 2. Israelachvili, J. N., Marcelja, S. & Horn, R. G. (1980) Q. Rev. % Biophys. 13, 121-200. % 3. Israelachvili, J. N. (1985) Intermolecular and Surface Forces % (Academic, New York), pp. 229-271. % \DLA{This section needs to be added still.} % PE = 1.33 (Kumar91) % CHOL = 1.21 (Kumar91) % PC=0.8 (Kumar91) % SM=0.8 (assumed by rz same as PC) % PS=1 (no refs so far; should be close to unity; rz) % % From Epand94LIB: a considerable amount of heat is evolved per mole of % bilayer stabilizer when such molecules are inserted into membranes % which are under a high curvature strain. If this energy were coupled % to a membrane event, such as the conformational change in a protein, % it could be the driving force responsible for such processes. % % From Epand94LIB: considerable energy may be released upon the % incorporation of certain molecules into membranes which have a low % radius of spontaneous curvature. {discuss with DL in terms of % catalytic activity; rz; 8/24/2010}. % % Kleinfeld93 showed that the transfer parameters of long-chain FFA % between the lipid vesicles depend on vesicle curvature and % composition. Transfer of stearic acid is much slower from LUV as % compared to SUV. % % Substitution of DPPC SUV for POPC SUV as the donor matrix revealed % some differences in the energetics of transfer that are associated % with the physical state of the donor particle. Whereas the free % energies of transfer of oleic acid from DPPC SUV and POPC SUV are % similar, the free energy of transfer of stearic acid from DPPC SUV is % much higher than that observed from POPC. % % From Kumar91: S is given by S = V/al, where V is the hydrocarbon % volume, a is the area of head group, and l is the critical length of % the hydrocarbon chain (1-3). a, V, and I are all estimable or % measurable (4), and the value of S can be calculated. The value of S % determines the aggregate formed by lipids or any amphiphiles upon % hydration. It has been shown that lipids aggregate to form spherical % micelles (S < 1/3), nonspherical (cylindrical) micelles (1/3 < S < % 1/2), bilayers (1/2 < S < 1), and reverse micelles or hexagonal (HII,) % phases (S > 1). However, theoreticians caution that the above % predicted limits set on the values of S are relatively insensitive to % the exact values of V and a but are strongly dependent upon the choice % of l (5). % % From Kumar91: Cylindrical and wedge-shaped molecules have been shown % to be essential for spontaneous vesiculation (6) as well as bilayer % stabilization (7, 8). % % From Kumar91: It has also been shown experimentally that molecules % having complementary molecular shapes can form bilayer structures. % Inverted cone (wedge)- shaped lysoPC and cone-shaped Chol or % unsaturated PE can interact stoichiometrically to form bilayer % structures (9-12). {quote to justify our formula for complementary % cones; rz} \section{Kinetic adjustments} \label{sec:kinetic_adjustments} In the mass action equation used in the formalism, both the forward and backwards kinetic constants ($k_\mathrm{f}$ and $k_\mathrm{b}$) are adjusted (by $k_{\mathrm{f}i\mathrm{adj}}$ and $k_{\mathrm{b}i\mathrm{adj}}$) to reflect the properties of the vesicle and the properties of the monomer entering or exiting the vesicle. \subsection{Forward Kinetic Adjustments} The forward rate constant adjustment, $k_{\mathrm{f}i\mathrm{adj}}$ takes into account unsaturation ($un_\mathrm{f}$), charge ($ch_\mathrm{f}$), curvature ($cu_\mathrm{f}$), length ($l_\mathrm{f}$), and complex formation ($CF1_\mathrm{f}$), each of which is modified depending on the specific component and the vesicle into which the component is entering; the final $k_{\mathrm{f}i\mathrm{adj}}$ is the product of the adjustments made separately for each property. \begin{equation} 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} \label{eq:kf_adj} \end{equation} \subsubsection{Unsaturation Forward} In order for a lipid to be inserted into a membrane, a void has to be formed for it to fill. Voids can be generated by the combination of unsaturated and saturated lipids forming heterogeneous domains. Void formation is increased when the unsaturation of lipids in the vesicle is widely distributed; in other words, the insertion of lipids into the membrane is greater when the standard deviation of the unsaturation is larger % %%% \RZ{May need to site (at least for us) works showing %%% mismatch-dependent ``defects''}. % Assuming that an increase in width of the distribution linearly decreases the free energy of activation, the $un_\mathrm{f}$ parameter must follow $a^{\mathrm{stdev}\left(un_\mathrm{vesicle}\right)}$ where $a > 1$, so a convenient starting base for $a$ which leads to a reasonable change in Eyring activation energy is $2$: \begin{equation} un_\mathrm{f} = 2^{\mathrm{stdev}\left(un_\mathrm{vesicle}\right)} \label{eq:unsaturation_forward} \end{equation} The mean $\mathrm{stdev}\left(un_\mathrm{vesicle}\right)$ in our simulations is around $1.5$, which leads to a $\Delta \Delta G^\ddagger$ of $\Sexpr{to.latex(format(digits=3,to.kcal(2^1.5)))} \frac{\mathrm{kcal}}{\mathrm{mol}}$, and a total range of possible values depicted in \cref{fig:unf_graph}. % \RZ{Explain here, or even earlier that the formulas were ad hoc % adjusted to correspond to ``reasonable'' changes in the Eyring % activation energies.} % It is not clear that the unsaturation of the inserted monomer will % affect the rate of the insertion positively or negatively, so we do % not include a term for it in this formalism. \begin{figure} <>= layout(matrix(1:2,nrow=1,ncol=2)) curve(2^x,from=0,to=sd(c(0,4)), xlab=expression(paste(stdev,group("(",italic(un[vesicle]),")"))), ylab=expression(italic(un[f]))) vps <- baseViewports() pushViewport(vps$figure) grid.text("A", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport() curve(to.kcal(2^x),from=0,to=sd(c(0,4)), xlab=expression(paste(stdev,group("(",italic(un[vesicle]),")"))), ylab="Activation Energy (kcal/mol)") vps <- baseViewports() pushViewport(vps$figure) grid.text("B", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport() @ \caption{Change in unsaturation forward ($un_\mathrm{f}$) (A) or activation energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the standard deviation of the unsaturation of the vesicle ($\mathrm{stdev}\left(un_\mathrm{vesicle}\right)$).} \label{fig:unf_graph} \end{figure} \subsubsection{Charge Forward} A charged lipid such as PS approaching a vesicle with an average charge of the same sign will experience repulsion, whereas those with different signs will experience attraction, the degree of which is dependent upon the charge of the monomer and the average charge of the vesicle. If either the vesicle or the monomer has no charge, there should be no effect of charge upon the rate. This leads us to the following equation, $a^{-\left ch_m}$ (which is similar to the numerator of Coulomb's law), where $\left$ is the average charge of the vesicle, and $ch_m$ is the charge of the monomer. If either $\left$ or $ch_m$ is 0, the adjustment parameter is 1 (no change), whereas it decreases if both are positive or negative, as the product of two real numbers with the same sign is always positive. The base for $a$ ($60$) was chosen ad hoc to correspond to a maximum of about $0.5\frac{\mathrm{kcal}}{\mathrm{mol}}$ change in activation energy for the common case, resulting in the following equation: \begin{equation} ch_\mathrm{f} = 60^{-\left<{ch}_v\right> {ch}_m} \label{eq:charge_forward} \end{equation} The most common $\left<{ch}_v\right>$ is around $-0.165$, which leads to a range of $\Delta \Delta G^\ddagger$ from $\Sexpr{format(digits=3,to.kcal(60^(-.165*-1)))} \frac{\mathrm{kcal}}{\mathrm{mol}}$ to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$, and the total range of possible values seen in \cref{fig:chf_graph}. \begin{figure} <>= trellis.device(new=F,color=TRUE) trellis.par.set(list(axis.line =list(col="transparent"))) pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off")) pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off")) x <- seq(-1,0,length.out=20) y <- seq(-1,0,length.out=20) grid <- expand.grid(x=x,y=y) grid$z <- as.vector(60^(-outer(x,y))) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, box=NA, scales=list(arrows=F,col=1,distance=0.75), xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30), ylab=list(label=expression(italic(ch[monomer])),rot=-35), zlab=list(label=expression(italic(ch[f])),rot=93)),newpage=FALSE) rm(x,y,grid) grid.text("A", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport(1) pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off")) x <- seq(-1,0,length.out=20) y <- seq(-1,0,length.out=20) grid <- expand.grid(x=x,y=y) grid$z <- as.vector(to.kcal(60^(-outer(x,y)))) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, strip=FALSE, scales=list(arrows=F,col=1,distance=0.75), xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30), ylab=list(label=expression(italic(ch[monomer])),rot=-35), zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE) rm(x,y,grid) grid.text("B", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport(2) @ \caption{Change in charge forward ($ch_\mathrm{f}$) (A) or activation energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the charge of the monomer entering ($ch_\mathrm{monomer}$) and the average charge of the vesicle ($\left$).} \label{fig:chf_graph} \end{figure} \subsubsection{Curvature Forward} Curvature is a measure of the intrinsic propensity of specific lipids to form micelles (positive curvature), inverted micelles (negative curvature), or planar sheets (neutral curvature)~\citep{Israelachvili1975:amphiphile_self_assembly}. % In this formalism, curvature is measured as the ratio of the volume of the lipid to the area of the head times the length of the lipid ($S=\frac{v}{l_c a}$), so negative curvature is bounded by $(0,1)$, neutral curvature is 1, and positive curvature is bounded by $(1,\infty)$. The curvature can be transformed using $\log$, which has the property of making the range of positive and negative curvature equal, and distributed about 0. As in the case of unsaturation, void formation is increased by the presence of lipids with mismatched curvature. Thus, a larger distribution of curvature in the vesicle increases the rate of lipid insertion into the vesicle. However, a component with curvature $cu^{-1}$ will cancel out a component with curvature $cu$, so we have to log transform (turning these into $-\log cu$ and $\log cu$), then take the absolute value ($\log cu$ and $\log cu$), and finally measure the width of the distribution, which in the case of exactly mismatched curvatures is $0$. Thus, by using the log transform to make the range of the lipid curvature equal between positive and negative, and taking the average to cancel out exactly mismatched curvatures, we come to an equation with the shape $a^{\left<\log cu_\mathrm{vesicle}\right>}$. An ad hoc base for $a$ which for the common case is $10$, yielding: \begin{equation} % cu_\mathrm{f} = 10^{\mathrm{stdev}\left|\log cu_\mathrm{vesicle}\right|} cu_\mathrm{f} = 10^{\left|\left<\log cu_\mathrm{monomer} \right>\right|\mathrm{stdev} \left|\log cu_\mathrm{vesicle}\right|} \label{eq:curvature_forward} \end{equation} The most common $\left|\left<\log {cu}_v\right>\right|$ is around $0.013$, which with the most common $\mathrm{stdev} \log cu_\mathrm{vesicle}$ of $0.213$ leads to a $\Delta \Delta G^\ddagger$ of $\Sexpr{format(digits=3,to.kcal(10^(0.13*0.213)))} \frac{\mathrm{kcal}}{\mathrm{mol}}$. This is a consequence of the relatively matched curvatures in our environment. The full range of $cu_\mathrm{f}$ values possible are shown in \cref{fig:cuf_graph}. % 1.5 to 0.75 3 to 0.33 \begin{figure} <>= trellis.device(new=F,color=TRUE) trellis.par.set(list(axis.line =list(col="transparent"))) pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off")) pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off")) grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))), sd(abs(log(c(1,0.33)))), sd(abs(log(c(0.33,3)))))),length.out=20), y=seq(0,max(c(mean(log(c(1,3)), mean(log(c(1,0.33))), mean(log(c(0.33,3)))))),length.out=20)) grid$z <- 10^(grid$x*grid$y) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, box=NA, scales=list(arrows=F,col=1,distance=0.75), xlab=list(label=expression(paste(stdev,group("(",italic(log(cu[vesicle])),")"))),rot=30), ylab=list(label=expression(paste(symbol("\341"),italic(log(cu[vesicle])),symbol("\361"))),rot=-35), zlab=list(label=expression(italic(cu[f])),rot=93)),newpage=FALSE) rm(grid) grid.text("A", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport(1) pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off")) grid <- expand.grid(x=seq(0,max(c(sd(abs(log(c(1,3)))), sd(abs(log(c(1,0.33)))),sd(abs(log(c(0.33,3)))))),length.out=20), y=seq(0,max(c(mean(log(c(1,3)), mean(log(c(1,0.33))), mean(log(c(0.33,3)))))),length.out=20)) grid$z <- to.kcal(10^(grid$x*grid$y)) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, box=NA, scales=list(arrows=F,col=1,distance=0.75), xlab=list(label=expression(paste(stdev,group("(",italic(log(cu[vesicle])),")"))),rot=30), ylab=list(label=expression(paste(symbol("\341"),italic(log(cu[vesicle])),symbol("\361"))),rot=-35), zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE) rm(grid) grid.text("B", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport(2) @ \caption{Change in curvature forward ($cu_\mathrm{f}$) (A) or activation energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the standard deviation of the log of the curvature of the vesicle ($\mathrm{stdev}\left(\log cu_\mathrm{vesicle}\right)$) and the mean of the log of the curvature of the vesicle ($\left<\log cu_\mathrm{vesicle}\right>$).} \label{fig:cuf_graph} \end{figure} \subsubsection{Length Forward} As in the case of unsaturation, void formation is easier when vesicles are made up of components of mismatched lengths. Thus, when the width of the distribution of lengths is larger, the forward rate should be greater as well, leading us to an equation of the form $x^{\mathrm{stdev} l_\mathrm{vesicle}}$, where $\mathrm{stdev} l_\mathrm{vesicle}$ is the standard deviation of the length of the components of the vesicle, which has a maximum possible value of $\Sexpr{format(digits=3,sd(c(rep(12,50),rep(24,50))))}$ and a minimum of 0 in this set of simulations. Based on activation energy, a reasonable base for $x$ is 2, leading to: \begin{equation} l_\mathrm{f} = 2^{\mathrm{stdev} l_\mathrm{vesicle}} \label{eq:length_forward} \end{equation} The most common $\mathrm{stdev} l_\mathrm{vesicle}$ is around $3.4$, which leads to a range of $\Delta \Delta G^\ddagger$ of $\Sexpr{format(digits=3,to.kcal(2^(3.4)))} \frac{\mathrm{kcal}}{\mathrm{mol}}$. % While it could be argued that increased length of the monomer could % affect the rate of insertion into the membrane, it's not clear whether % it would increase (by decreasing the number of available hydrogen % bonds, for example) or decrease (by increasing the time taken to fully % insert the acyl chain, for example) the rate of insertion or to what % degree, so we do not take it into account in this formalism. % \DLA{Incorporate McLean84 here} % From McLean84LIB: Although it is difficult to measure cmc values for % the sparingly soluble lipids used in this study, estimates for % lysopalmitoylphosphatidylcholine( 7 X l0-6 M; Haberland \& Reynolds, % 1975), cholesterol (12.1 X 10-8 M, extrapolated to infinite dilution; % Haberland \& Reynolds, 1973), and dipalmitoylphosphatidylcholine (4.6 % X l0-10 M; Smith \& Tanford, 1972) are available. A value of 1.1 X % 10-8 M for DMPC was estimated from the linear relationship between ln % cmc and the number of carbons in the PC acyl chain by using data for n % = 7, 8, 10, and 16 [summarized in Tanford (1980)]. % % From Nichols85: The magnitude of the dissociation rate constant % decreases by a factor of approximately 3.2 per carbon increase in acyl % chain length (see Table II here) {take into account for the formula; % rz 8/17/2010}. % % From Nichols85: The magnitude of the partition coefficient increases % with acyl chain length [Keq(M-C6-NBD-PC) = (9.8±2.1) X 106 M-1 and Keq % (P-C6-NBD-PC) = (9.4±1.3) x 107 M-1. {take into account for the % formula; rz 8/17/2010}. % % From Nichols85: The association rate constant is independent of acyl % chain length. {take into account for the formula; rz 8/17/2010}. \begin{figure} <>= layout(matrix(1:2,nrow=1,ncol=2)) curve(2^x,from=0,to=sd(c(12,24)), xlab=expression(paste(stdev,group("(",italic(l[vesicle]),")"))), ylab=expression(italic(l[f])), las=1 ) vps <- baseViewports() pushViewport(vps$figure) grid.text("A", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport() curve(to.kcal(2^x),from=0,to=sd(c(12,24)), xlab=expression(paste(stdev,group("(",italic(l[vesicle]),")"))), ylab="Activation Energy (kcal/mol)", las=1) vps <- baseViewports() pushViewport(vps$figure) grid.text("B", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport() @ \caption{Change in length forward ($l_\mathrm{f}$) (A) or activation energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the standard deviation of the length of the vesicle ($\mathrm{stdev}\left(l_\mathrm{vesicle}\right)$).} \label{fig:lf_graph} \end{figure} \subsubsection{Complex Formation} There is no contribution of complex formation to the forward reaction rate in the current formalism. \begin{equation} CF1_\mathrm{f}=1 \label{eq:complex_formation_forward} \end{equation} \subsection{Backward adjustments ($k_{\mathrm{b}i\mathrm{adj}}$)} Just as the forward rate constant adjustment $k_{\mathrm{f}i\mathrm{adj}}$ does, the backwards rate constant adjustment $k_{\mathrm{b}i\mathrm{adj}}$ takes into account unsaturation ($un_\mathrm{b}$), charge ($ch_\mathrm{b}$), curvature ($cu_\mathrm{b}$), length ($l_\mathrm{b}$), and complex formation ($CF1_\mathrm{b}$), each of which are modified depending on the specific component and the vesicle from which the component is exiting: \begin{equation} 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} \label{eq:kb_adj} \end{equation} \subsubsection{Unsaturation Backward} Unsaturation also influences the ability of a lipid molecule to leave a membrane. If a molecule has an unsaturation level which is different from the surrounding membrane, it will be more likely to leave the membrane. The more different the unsaturation level is, the greater the propensity for the lipid molecule to leave. However, a vesicle with some unsaturation (eg. average unsaturation of 2) is more favorable for lipids with more unsaturation (eg. 3) than lipids with equivalently less unsaturation (eg. 1), so the difference in energy between unsaturation is not linear. Therefore, an equation with the shape $x^{\left| y^{-\left< un_\mathrm{vesicle}\right> }-y^{-un_\mathrm{monomer}} \right| }$, where $\left$ is the average unsaturation of the vesicle and $un_\mathrm{monomer}$ is the unsaturation of the monomer, allows for increasing the efflux of molecules from membranes where they strongly mismatch, while allowing vesicles with greater unsaturation to tolerate greater unsaturation mismatch between monomers and the vesicle. The bases x and y were chosen ad hoc to produce reasonable Eyring energies of activation over the range of unsaturations expected, leading to: \begin{equation} un_\mathrm{b} = 7^{1-\left(20\left(2^{-\left} - 2^{-un_\mathrm{monomer}} \right)^2+1\right)^{-1}} \label{eq:unsaturation_backward} \end{equation} The most common $\left$ is around $1.7$, which leads to a range of $\Delta \Delta G^\ddagger$ from $\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-0)^2+1))))} \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with 0 unsaturation to $\Sexpr{format(digits=3,to.kcal(7^(1-1/(5*(2^-1.7-2^-4)^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with 4 unsaturations. See \cref{fig:unb_graph} for the full range of possible values. \begin{figure} <>= trellis.device(new=F,color=TRUE) trellis.par.set(list(axis.line =list(col="transparent"))) pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off")) pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off")) grid <- expand.grid(x=seq(0,4,length.out=20), y=seq(0,4,length.out=20)) grid$z <- (7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1))) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, box=NA, scales=list(arrows=F,col=1,distance=0.75), xlab=list(label=expression(paste(symbol("\341"),italic(un[vesicle]),symbol("\361"))),rot=30), ylab=list(label=expression(italic(un[monomer])),rot=-35), zlab=list(label=expression(italic(un[b])),rot=93)),newpage=FALSE) grid.text("A", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport(1) pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off")) grid <- expand.grid(x=seq(0,4,length.out=20), y=seq(0,4,length.out=20)) grid$z <- to.kcal((7^(1-1/(5*(2^-grid$x-2^-grid$y)^2+1)))) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, strip=FALSE, scales=list(arrows=F,col=1,distance=0.75), xlab=list(label=expression(paste(symbol("\341"),italic(un[vesicle]),symbol("\361"))),rot=30), ylab=list(label=expression(italic(un[monomer])),rot=-35), zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE) rm(grid) grid.text("B", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport(2) @ \caption{Change in unsaturation backward ($un_\mathrm{b}$) (A) or activation energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the unsaturation of the monomer leaving ($un_\mathrm{monomer}$) and the average unsaturation of the vesicle ($\left$).} \label{fig:unb_graph} \end{figure} \subsubsection{Charge Backwards} As in the case of monomers entering a vesicle, monomers leaving a vesicle leave faster if their charge has the same sign as the average charge vesicle. An equation of the form $ch_\mathrm{b} = a^{\left ch_m}$ is then appropriate, and using a base of $a=20$ yields: \begin{equation} ch_\mathrm{b} = 20^{\left<{ch}_v\right> {ch}_m} \label{eq:charge_backwards} \end{equation} The most common $\left$ is around $-0.164$, which leads to a range of $\Delta \Delta G^\ddagger$ from $\Sexpr{format(digits=3,to.kcal(20^(-.164*-1)))} \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with charge $-1$ to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with charge $0$. See \cref{fig:chb_graph} for the full range of possible values of $ch_\mathrm{b}$. \begin{figure} <>= trellis.device(new=F,color=TRUE) trellis.par.set(list(axis.line =list(col="transparent"))) pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off")) pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off")) x <- seq(-1,0,length.out=20) y <- seq(-1,0,length.out=20) grid <- expand.grid(x=x,y=y) grid$z <- as.vector(20^(outer(x,y))) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, box=NA, scales=list(arrows=F,col=1,distance=0.75), xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30), ylab=list(label=expression(italic(ch[monomer])),rot=-35), zlab=list(label=expression(italic(ch[b])),rot=93)),newpage=FALSE) rm(x,y,grid) grid.text("A", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport(1) pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off")) x <- seq(-1,0,length.out=20) y <- seq(-1,0,length.out=20) grid <- expand.grid(x=x,y=y) grid$z <- to.kcal(as.vector(20^(outer(x,y)))) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, strip=FALSE, scales=list(arrows=F,col=1,distance=0.75), xlab=list(label=expression(paste(symbol("\341"),italic(ch[vesicle]),symbol("\361"))),rot=30), ylab=list(label=expression(italic(ch[monomer])),rot=-35), zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE) rm(x,y,grid) grid.text("B", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport(2) @ \caption{Change in charge backward ($ch_\mathrm{b}$) (A) or activation energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the charge of the monomer leaving ($ch_\mathrm{monomer}$) and the average charge of the vesicle ($\left$).} \label{fig:chb_graph} \end{figure} \subsubsection{Curvature Backwards} The less a monomer's intrinsic curvature matches the average curvature of the vesicle in which it is in, the greater its rate of efflux. If the difference is 0, $cu_\mathrm{f}$ needs to be one. To map negative and positive curvature to the same range, we also need take the logarithm. Positive and negative curvature lipids cancel each other out, resulting in an average curvature of 0; the average of the log also has this property. Increasing mismatches in curvature increase the rate of efflux, but asymptotically. An equation which satisfies these criteria has the form $cu_\mathrm{f} = a^{1-\left(b\left( \left< \log cu_\mathrm{vesicle} \right> -\log cu_\mathrm{monomer}\right)^2+1\right)^{-1}}$. An alternative form would use the absolute value of the difference, however, this yields a cusp and sharp increase about the curvature equilibrium. We have chosen bases of $a=7$ and $b=20$ ad hoc, based on activation energy considerations. \begin{equation} cu_\mathrm{b} = 7^{1-\left(20\left(\left< \log cu_\mathrm{vesicle} \right> -\log cu_\mathrm{monomer} \right)^2+1\right)^{-1}} \label{eq:curvature_backwards} \end{equation} The most common $\left<\log cu_\mathrm{vesicle}\right>$ is around $-0.013$, which leads to a range of $\Delta \Delta G^\ddagger$ from $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(0.8))^2+1))))} \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature 0.8 to $\Sexpr{format(digits=3,to.kcal(7^(1-1/(20*(-0.013-log(1.3))^2+1))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature 1.3 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with curvature near 1. The full range of values possible for $cu_\mathrm{b}$ are shown in \cref{fig:cub_graph} % \RZ{What about the opposite curvatures that actually do fit to each % other?} % From Kumar91: It has also been shown experimentally that molecules % having complementary molecular shapes can form bilayer structures. % Inverted cone (wedge)- shaped lysoPC and cone-shaped Chol or % unsaturated PE can interact stoichiometrically to form bilayer % structures (9-12). {quote to justify our formula for complementary % cones; rz} \begin{figure} <>= trellis.device(new=F,color=TRUE) trellis.par.set(list(axis.line =list(col="transparent"))) pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off")) pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off")) grid <- expand.grid(x=log(seq(0.8,1.33,length.out=20)), y=seq(0.8,1.33,length.out=20)) grid$z <- 7^(1-1/(20*(grid$x-log(grid$y))^2+1)) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, box=NA, scales=list(arrows=F,col=1,distance=0.75), xlab=list(label=expression(paste(symbol("\341"),log(italic(cu[vesicle])),symbol("\361"))),rot=30), ylab=list(label=expression(italic(cu[monomer])),rot=-35), zlab=list(label=expression(italic(cu[b])),rot=93)),newpage=FALSE) rm(grid) grid.text("A", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport(1) pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off")) grid <- expand.grid(x=log(seq(0.8,1.33,length.out=20)), y=seq(0.8,1.33,length.out=20)) grid$z <- to.kcal(7^(1-1/(20*(grid$x-log(grid$y))^2+1))) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, strip=FALSE, scales=list(arrows=F,col=1,distance=0.75), xlab=list(label=expression(paste(symbol("\341"),log(italic(cu[vesicle])),symbol("\361"))),rot=30), ylab=list(label=expression(italic(cu[monomer])),rot=-35), zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE) rm(grid) grid.text("B", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport(2) @ \caption{Change in curvature backward ($cu_\mathrm{b}$) (A) or activation energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the curvature of the monomer leaving ($ch_\mathrm{monomer}$) and the average of the log of the curvature of the vesicle ($\left<\log cu_\mathrm{vesicle}\right>$).} \label{fig:cub_graph} \end{figure} \subsubsection{Length Backwards} In a model membrane, the dissociation constant increases by a factor of approximately 3.2 per carbon decrease in acyl chain length~% \citep{Nichols1985:phospholipid_monomer_vesicle_thermodynamics}. Unfortunately, the experimental data known to us only measures the effect of chain lengths less than or equal to the bulk lipid, not the effect of lengths exceeding it, and is only known for one bulk lipid component (DOPC). We assume then, that the increase is in relationship to the average vesicle, and that lipids with larger acyl chain length will also show an increase in their dissociation constant. \begin{equation} l_\mathrm{b} = 3.2^{\left|\left-l_\mathrm{monomer}\right|} \label{eq:length_backward} \end{equation} The most common $\left$ is around $17.75$, which leads to a range of $\Delta \Delta G^\ddagger$ from $\Sexpr{format(digits=3,to.kcal(3.2^abs(12-17.75)))} \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with length 12 to $\Sexpr{format(digits=3,to.kcal(3.2^abs(24-17.75)))}\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with length 24 to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with length near 18. The full range of possible values of $l_\mathrm{b}$ are shown in \cref{fig:lb_graph} % (for methods? From McLean84LIB: The activation free energies and free % energies of transfer from self-micelles to water increase by 2.2 and % 2.1 kJ mol-' per methylene group, respectively. \RZ{see if we can use it % to justify arranging our changed in activating energy around 1 % kcal/mol}). \begin{figure} <>= trellis.device(new=F,color=TRUE) trellis.par.set(list(axis.line =list(col="transparent"))) pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off")) pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off")) grid <- expand.grid(x=seq(12,24,length.out=20), y=seq(12,24,length.out=20)) grid$z <- 3.2^(abs(grid$x-grid$y))/1e6 print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, box=NA, scales=list(arrows=F,col=1,distance=0.75), xlab=list(label=expression(paste(symbol("\341"),italic(l[vesicle]),symbol("\361"))),rot=30), ylab=list(label=expression(italic(l[monomer])),rot=-35), zlab=list(label=expression(italic(l[b])%*%10^-6),rot=93)),newpage=FALSE) rm(grid) grid.text("A", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport(1) pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off")) grid <- expand.grid(x=seq(12,24,length.out=20), y=seq(12,24,length.out=20)) grid$z <- to.kcal(3.2^(abs(grid$x-grid$y))) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, strip=FALSE, scales=list(arrows=F,col=1,distance=0.75), xlab=list(label=expression(paste(symbol("\341"),italic(l[vesicle]),symbol("\361"))),rot=30), ylab=list(label=expression(italic(l[monomer])),rot=-35), zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE) rm(grid) grid.text("B", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport(2) @ \caption{Change in length backwards ($l_\mathrm{b}\times 10^{-6}$) (A) or activation energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the length of the monomer leaving ($l_\mathrm{monomer}$) and the average length of lipids in the vesicle ($\left$).} \label{fig:lb_graph} \end{figure} \subsubsection{Complex Formation Backward} Complex formation ($CF1$) describes the interaction between CHOL and PC or SM, where PC or SM protects the hydroxyl group of CHOL from interactions with water % \citep{Huang1999:chol_solubility_model, Huang1999:cholesterol_solubility, McConnell2006, McConnell2003}% . PC ($CF1=2$) can interact with two CHOL, and SM ($CF1=3$) with three CHOL ($CF1=-1$). If the average of $CF1$ is positive (excess of SM and PC with regards to complex formation), components with negative $CF1$ (CHOL) will be retained. If average $CF1$ is negative, components with positive $CF1$ are retained. An equation which has this property is $CF1_\mathrm{b}=a^{\left CF1_\mathrm{monomer}-\left|\left CF1_\mathrm{monomer}\right|}$, where difference of the exponent is zero if the average $CF1$ and the $CF1$ of the component have the same sign, or double the product if the signs are different. Based on activation energy considerations, we took an ad hoc base for $a$ as $1.5$. \begin{equation} CF1_\mathrm{b}=1.5^{\left CF1_\mathrm{monomer}-\left|\left CF1_\mathrm{monomer}\right|} \label{eq:complex_formation_backward} \end{equation} The most common $\left$ is around $0.925$, which leads to a range of $\Delta \Delta G^\ddagger$ from $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*-1-abs(0.925*-1))))} \frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex formation $-1$ to $\Sexpr{format(digits=3,to.kcal(1.5^(0.925*2-abs(0.925*2))))}\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex formation $2$ to $0\frac{\mathrm{kcal}}{\mathrm{mol}}$ for monomers with complex formation $0$. The full range of possible values for $CF1_\mathrm{b}$ are depicted in \cref{fig:cf1b_graph}. \begin{figure} <>= trellis.device(new=F,color=TRUE) trellis.par.set(list(axis.line =list(col="transparent"))) pushViewport(viewport(layout=grid.layout(nrow=1,ncol=2),clip="off")) pushViewport(viewport(layout.pos.row=1,layout.pos.col=1,clip="off")) grid <- expand.grid(x=seq(-1,3,length.out=20), y=seq(-1,3,length.out=20)) grid$z <- 1.5^(grid$x*grid$y-abs(grid$x*grid$y)) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, box=NA, scales=list(arrows=F,col=1,distance=0.75), xlab=list(label=expression(paste(symbol("\341"),italic(CF1[vesicle]),symbol("\361"))),rot=30), ylab=list(label=expression(italic(CF1[monomer])),rot=-35), zlab=list(label=expression(italic(CF1[b])),rot=93)),newpage=FALSE) rm(grid) grid.text("A", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport(1) pushViewport(viewport(layout.pos.row=1,layout.pos.col=2,clip="off")) grid <- expand.grid(x=seq(-1,3,length.out=20), y=seq(-1,3,length.out=20)) grid$z <- to.kcal(1.5^(grid$x*grid$y-abs(grid$x*grid$y))) print(wireframe(z~x*y,grid,cuts=50, drape=TRUE, strip=FALSE, scales=list(arrows=F,col=1,distance=0.75), xlab=list(label=expression(paste(symbol("\341"),italic(CF1[vesicle]),symbol("\361"))),rot=30), ylab=list(label=expression(italic(CF1[monomer])),rot=-35), zlab=list(label="Activation Energy (kcal/mol)",rot=93)),newpage=FALSE) rm(grid) grid.text("B", just=c("left","top"), gp=gpar(fontsize=48), y=unit(1,"npc")-unit(0.25,"lines"), x=unit(0,"npc")+unit(0.25,"lines")) popViewport(2) @ \caption{Change in complex formation 1 backwards ($CF1_\mathrm{b}$) (A) or activation energy in $\frac{\mathrm{kcal}}{\mathrm{mol}}$ (B) versus the complex formation of the monomer leaving ($CF1_\mathrm{monomer}$) and the average complex formation of the vesicle ($\left$).} \label{fig:cf1b_graph} \end{figure} \section{Simulation Methodology} \subsection{Overall Architecture} The simulations are currently run using a program written in perl with various modules to handle the subsidiary parts. It produces output for each generation, including the step immediately preceding and immediately following a vesicle split (and optionally, each step) that is written to a state file which contains the state of the vesicle, environment, kinetic parameters, program invocation options, time, and various other parameters necessary to recreate the state vector at a given time. This output file is then read by a separate program in perl to produce different output which is generated out-of-band; output which includes graphs and statistical analysis is performed using R~\citep{R:maincite} (and various grid graphics modules) called from the perl program. The separation of simulation and output generation allows refining output, and simulations performed on different versions of the underlying code to be compared using the same analysis methods and code. It also allows later simulations to be restarted from a specific generation, utilizing the same environment. Random variables of different distributions are calculated using Math::Random (0.71), which is seeded for each run using entropy from the linux kernel's urandom device. Code is available upon request. \subsection{Environment Creation} \subsubsection{Components} The environment contains concentrations of components. In the current set of simulations, there are \Sexpr{1+4*5*7} different possible components, consisting of PC, PE, PS, SM, and CHOL, with all lipids except for CHOL having 5 possible unsaturations ranging from 0 to 4, and 7 lengths from $12,14,...,22$ ($7\cdot 5\cdot4+1=\Sexpr{1+4*5*7}$). In cases where the environment has less than the maximum number of components, the components are selected in order without replacement from a randomly shuffled deck of components (with the exception of CHOL) represented once until the desired number of unique components is obtained. CHOL is overrepresented $\Sexpr{5*7}$ times to be at the level of other lipid types, ensuring that the probability of CHOL being absent in the environment is the same as the probability of any one of the other lipid types (PS, PE, etc.) being absent. This reduces the number of simulations with a small number of components which are devoid of CHOL. \subsubsection{Concentrations} Once the components of the environment have been selected, their concentrations are determined. In experiments where the environmental concentration is the same across all lipid components, the concentration is set to $10^{-10}\mathrm{M}$. In other cases, the environmental concentration is set to a random number from a gamma distribution with shape parameter 2 and an average of $10^{-10}\mathrm{M}$. The base concentration ($10^{-10}\mathrm{M}$) can also be altered at the initialization of the experiment to specific values for specific lipid types or components. \subsection{Initial Vesicle Creation} Initial vesicles are seeded by selecting lipid molecules from the environment until the vesicle reaches a specific starting size. The vesicle starting size has gamma distribution with shape parameter 2 and a mean of the per-simulation specified starting size, with a minimum of 5 lipid molecules. Lipid molecules are then selected to be added to the lipid membrane according to four different methods. In the constant method, molecules are added in direct proportion to their concentration in the environment. The uniform method adds molecules in proportion to their concentration in the environment scaled by a uniform random value, whereas the random method adds molecules in proportion to a uniform random value. The final method is a binomial method, which adjust the probability of adding a molecule of a specific component by the concentration of that component, and then adds the molecules one by one to the membrane. This final method is also used in the first three methods to add any missing molecules to the starting vesicle which were unallocated due to the requirement to add integer numbers of molecules. (For example, if a vesicle was to have 10 molecules split evenly between three lipid types, the 10th molecule would be assigned by randomly choosing between the three lipid types, weighted by their concentration in the environment.) \subsection{Simulation Step} Once the environment has been created and the initial vesicle has been formed, molecules join and leave the vesicle based on the kinetic parameters and state equation discussed above until the vesicle splits forming two progeny vesicles, one of which the program continues to follow. \subsubsection{Calculation of Vesicle Properties} \label{sec:ves_prop} $S_\mathrm{vesicle}$ is the surface area of the vesicle, and is the sum of the surface area of all of the individual lipid components; each lipid type has a different surface area; we assume that the lipid packing is optimal, and there is no empty space. For the other vesicle properties (length, unsaturation, charge, and curvature), we calculate the mean, the standard deviation, the mean of the log, the mean of the absolute value of the log, the standard deviation of the log, and the standard deviation of the absolute value of the log. All cases, when we take the log, we take the log of the absolute value, and map $\log(0)$ to $0$. For the purpose of plotting vesicle properties, we only plot the mean of the property. \subsubsection{Joining and Leaving of Lipid Molecules} Determining the number of molecules to add to the lipid membrane ($n_i$) requires knowing $k_{\mathrm{f}i_\mathrm{adj}}$, the surface area of the vesicle $S_\mathrm{vesicle}$ (see \cref{sec:ves_prop}), the time interval $dt$ during which lipids are added, the base $k_{\mathrm{f}i}$, and the concentration of the monomer in the environment $\left[C_{i\mathrm{vesicle}}\right]$ (see \cref{eq:state}). $k_{\mathrm{f}i\mathrm{adj}}$ is calculated (see \cref{eq:kf_adj}) based on the vesicle properties and their hypothesized effect on the rate (in as many cases as possible, experimentally based) (see \cref{sec:kinetic_adjustments} for details). $dt$ can be varied (see \cref{sec:step_duration}), but for a given step is constant. This leads to the following: $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$ In the cases where $n_i > 1$, the integer number of molecules is added. Fractional $n_i$ or the fractional remainder after the addition of the integer molecules are the probability of adding a specific molecule, and are compared to a uniformly distributed random value between 0 and 1. If the random value is less than or equal to the fractional part of $n_i$, an additional molecule is added. Molecules leaving the vesicle are handled in a similar manner, with $n_i = k_{\mathrm{b}i}k_{\mathrm{b}i_\mathrm{adj}}C_{i_\mathrm{vesicle}}\mathrm{N_A}dt$. Before addition or removal of molecules, the properties of the vesicle are calculated; this is done so that molecules can be considered to be added and removed at the same instant, even though additions are handled first programmatically. This also avoids cases where a removal would have resulted in a negative number of molecules for a particular type. \subsubsection{Step duration} \label{sec:step_duration} $dt$ is the time taken for each step of joining, leaving, and checking split. In the current implementation, it starts with a value of $10^{-6}$~$\mathrm{s}$ but this is modified in between each step if the number of molecules joining or leaving is too large or small. If more than half of the starting vesicle size molecules join or leave in a single step, $dt$ is reduced by half. If less than the starting vesicle size molecules join or leave in 100 steps, $dt$ is doubled. This is necessary to curtail run times and to automatically adjust to different experimental runs. \subsubsection{Vesicle splitting} If a vesicle has grown to a size which is more than double the starting vesicle size, the vesicle splits. Molecules are assigned to the progeny vesicles at random, with each progeny vesicle having an equal chance of getting a single molecule. The number of molecules to assign to the first vesicle has binomial distribution with a probability of an event in each trial of 0.5 and a number of trials equal to the number of molecules. \subsection{Output} The environment, initial vesicle, and the state of the vesicle immediately before and immediately after splitting are stored to produce later output. % \section{Analyzing output} % % The analysis of output is handled by a separate perl program which % shares many common modules with the simulation program. Current output % includes simulation progress, summary tables, summary statistics, and % various graphs. % % \subsection{PCA plots} % % Two major groups of axes that contribute to vesicle variation between % subsequent generations are the components and properties of vesicles. % Each component in a vesicle is an axis on its own; it can be measured % either as an absolute number of molecules in each component, or the % fraction of molecules of that component over the total number of % molecules; the second approach is often more convenient, as it allows % vesicles with different number of molecules to be directly compared % (though it hides any effect of vesicle size). % % In order to visualize the variation between and transition of % subsequent generations of vesicles from their initial state in the % simulation, to their final state at the termination of the simulation, % we plot the projection of the generations onto the two principle PCA % axes (and alternatively, any pairing of the axes). % % \subsection{Carpet plots} % % Carpet plots show the distance/similarity between the composition or % properties of all generations in a single run. The difference between % each group of vesicle can be calculated using Euclidean distance % (properties and compositions) or H similarity (composition only). We % must use Euclidean distance for properties because the H distance % metric is invalid when comparing negative vector elements to positive % vector elements. % % In addition to showing distance or similarity, carpet plots also % depict propersomes and composomes as square boxes on the diagonals and % propertypes and compotypes as rectangles off the diagonals, each % propertype or compotype with a distinct color. % % \subsection{Propersomes, propertypes, composomes and compotypes} % % A generation is considered to be a propersome if it is less than $0.6$ % units (by Euclidean distance of normalized properties) away from the % generation immediately following and preceding. Likewise, a generation % is in a composome if its H similarity is more than $0.9$ (by the % normalized H metric) from the preceding and following generations. % Propersomes and composomes are then assigned to propertypes and % compotypes using Paritioning Around Medioids % (PAM). All values of $k$ between 2 and 15 % (or the number of propersomes and composomes, if that is less than 15) % are tried, and the clustering with the smallest % silhouette~\citep{Rousseeuw1987:silhouettes} is chosen as the ideal % clustering~\citep{Shenhav2005:pgard}. \bibliographystyle{unsrtnat} \bibliography{references.bib} \end{document}