ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mattDisertation/lipid.tex
(Generate patch)

Comparing trunk/mattDisertation/lipid.tex (file contents):
Revision 971 by mmeineke, Wed Jan 21 02:25:18 2004 UTC vs.
Revision 1092 by mmeineke, Tue Mar 16 21:35:16 2004 UTC

# Line 4 | Line 4
4  
5   \section{\label{lipidSec:Intro}Introduction}
6  
7 < \section{\label{lipidSec:Methods}Methods}
7 > In the past 10 years, increasing computer speeds have allowed for the
8 > atomistic simulation of phospholipid bilayers for increasingly
9 > relevant lengths of time.  These simulations have ranged from
10 > simulation of the gel ($L_{\beta}$) phase of
11 > dipalmitoylphosphatidylcholine (DPPC),\cite{lindahl00} to the
12 > spontaneous aggregation of DPPC molecules into fluid phase
13 > ($L_{\alpha}$) bilayers.\cite{marrink01} With the exception of a few
14 > ambitious
15 > simulations, \cite{marrink01:undulation,marrink:2002,lindahl00} most
16 > investigations are limited to a range of 64 to 256
17 > phospholipids.\cite{lindahl00,sum:2003,venable00,gomez:2003,smondyrev:1999,marrink01}
18 > The expense of the force calculations involved when performing these
19 > simulations limits the system size. To properly hydrate a bilayer, one
20 > typically needs 25 water molecules for every lipid, bringing the total
21 > number of atoms simulated to roughly 8,000 for a system of 64 DPPC
22 > molecules. Added to the difficulty is the electrostatic nature of the
23 > phospholipid head groups and water, requiring either the
24 > computationally expensive, direct Ewald sum or the slightly faster particle
25 > mesh Ewald (PME) sum.\cite{nina:2002,norberg:2000,patra:2003} These factors
26 > all limit the system size and time scales of bilayer simulations.
27  
28 < \subsection{\label{lipidSec:lipidMedel}The Lipid Model}
28 > Unfortunately, much of biological interest happens on time and length
29 > scales well beyond the range of current simulation technologies. One
30 > such example is the observance of a ripple ($P_{\beta^{\prime}}$)
31 > phase which appears between the $L_{\beta}$ and $L_{\alpha}$ phases of
32 > certain phospholipid bilayers
33 > (Fig.~\ref{lipidFig:phaseDiag}).\cite{katsaras00,sengupta00} These
34 > ripples are known from x-ray diffraction data to have periodicities on
35 > the order of 100-200~$\mbox{\AA}$.\cite{katsaras00} A simulation on
36 > this length scale would have approximately 1,300 lipid molecules with
37 > an additional 25 water molecules per lipid to fully solvate the
38 > bilayer. A simulation of this size is impractical with current
39 > atomistic models.
40  
41   \begin{figure}
42 + \centering
43 + \includegraphics[width=\linewidth]{ripple.eps}
44 + \caption[Diagram of the bilayer gel to fluid phase transition]{Diagram showing the $P_{\beta^{\prime}}$ phase as a transition between the $L_{\beta}$ and $L_{\alpha}$ phases}
45 + \label{lipidFig:phaseDiag}
46 + \end{figure}
47  
48 < \caption{Schematic diagram of the single chain phospholipid model}
48 > The time and length scale limitations are most striking in transport
49 > phenomena.  Due to the fluid-like properties of lipid membranes, not
50 > all small molecule diffusion across the membranes happens at pores.
51 > Some molecules of interest may incorporate themselves directly into
52 > the membrane.  Once there, they may exhibit appreciable waiting times
53 > (on the order of 10's to 100's of nanoseconds) within the
54 > bilayer. Such long simulation times are nearly impossible to obtain
55 > when integrating the system with atomistic detail.
56  
57 < \label{lipidFig:lipidModel}
57 > To address these issues, several schemes have been proposed.  One
58 > approach by Goetz and Lipowsky\cite{goetz98} is to model the entire
59 > system as Lennard-Jones spheres. Phospholipids are represented by
60 > chains of beads with the top most beads identified as the head
61 > atoms. Polar and non-polar interactions are mimicked through
62 > attractive and soft-repulsive potentials respectively.  A model
63 > proposed by Marrinck \emph{et. al}.\cite{marrink04}~uses a similar
64 > technique for modeling polar and non-polar interactions with
65 > Lennard-Jones spheres. However, they also include charges on the head
66 > group spheres to mimic the electrostatic interactions of the
67 > bilayer. The solvent spheres are kept charge-neutral and
68 > interact with the bilayer solely through an attractive Lennard-Jones
69 > potential.
70  
71 + The model used in this investigation adds more information to the
72 + interactions than the previous two models, while still balancing the
73 + need for simplification of atomistic detail.  The model uses
74 + unified-atom Lennard-Jones spheres for the head and tail groups of the
75 + phospholipids, allowing for the ability to scale the parameters to
76 + reflect various sized chain configurations while keeping the number of
77 + interactions small.  What sets this model apart, however, is the use
78 + of dipoles to represent the electrostatic nature of the
79 + phospholipids. The dipole electrostatic interaction is shorter range
80 + than Coulombic ($\frac{1}{r^3}$ versus $\frac{1}{r}$), and therefore
81 + eliminates the need for the costly Ewald sum.
82 +
83 + Another key feature of this model is the use of a dipolar water model
84 + to represent the solvent. The soft sticky dipole ({\sc ssd}) water
85 + \cite{liu96:new_model} relies on the dipole for long range electrostatic
86 + effects, but also contains a short range correction for hydrogen
87 + bonding. In this way the simulated systems in this research mimic the
88 + entropic contribution to the hydrophobic effect due to hydrogen-bond
89 + network deformation around a non-polar entity, \emph{i.e.}~the
90 + phospholipid molecules. This effect has been missing from previous
91 + reduced models.
92 +
93 + The following is an outline of this chapter.
94 + Sec.~\ref{lipidSec:Methods} is an introduction to the lipid model used
95 + in these simulations, as well as clarification about the water model
96 + and integration techniques. The various simulations explored in this
97 + research are outlined in
98 + Sec.~\ref{lipidSec:ExpSetup}. Sec.~\ref{lipidSec:resultsDis} gives a
99 + summary and interpretation of the results.  Finally, the conclusions
100 + of this chapter are presented in Sec.~\ref{lipidSec:Conclusion}.
101 +
102 + \section{\label{lipidSec:Methods}Methods}
103 +
104 + \subsection{\label{lipidSec:lipidModel}The Lipid Model}
105 +
106 + \begin{figure}
107 + \centering
108 + \includegraphics[width=\linewidth]{twoChainFig.eps}
109 + \caption[The two chained lipid model]{Schematic diagram of the double chain phospholipid model. The head group (in red) has a point dipole, $\boldsymbol{\mu}$, located at its center of mass. The two chains are eight methylene groups in length.}
110 + \label{lipidFig:lipidModel}
111   \end{figure}
112  
113   The phospholipid model used in these simulations is based on the
114   design illustrated in Fig.~\ref{lipidFig:lipidModel}.  The head group
115   of the phospholipid is replaced by a single Lennard-Jones sphere of
116 < diameter $fix$, with $fix$ scaling the well depth of its van der Walls
117 < interaction.  This sphere also contains a single dipole of magnitude
118 < $fix$, where $fix$ can be varied to mimic the charge separation of a
116 > diameter $\sigma_{\text{head}}$, with $\epsilon_{\text{head}}$ scaling
117 > the well depth of its van der Walls interaction.  This sphere also
118 > contains a single dipole of magnitude $|\boldsymbol{\mu}|$, where
119 > $|\boldsymbol{\mu}|$ can be varied to mimic the charge separation of a
120   given phospholipid head group.  The atoms of the tail region are
121 < modeled by unified atom beads.  They are free of partial charges or
122 < dipoles, containing only Lennard-Jones interaction sites at their
123 < centers of mass.  As with the head groups, their potentials can be
124 < scaled by $fix$ and $fix$.
121 > modeled by beads representing multiple methyl groups.  They are free
122 > of partial charges or dipoles, and contain only Lennard-Jones
123 > interaction sites at their centers of mass.  As with the head groups,
124 > their potentials can be scaled by $\sigma_{\text{tail}}$ and
125 > $\epsilon_{\text{tail}}$.
126  
127 < The long range interactions between lipids are given by the following:
127 > The possible long range interactions between atomic groups in the
128 > lipids are given by the following:
129   \begin{equation}
130 < EQ Here
130 > V_{\text{LJ}}(r_{ij}) =
131 >        4\epsilon_{ij} \biggl[
132 >        \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12}
133 >        - \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6}
134 >        \biggr]
135   \label{lipidEq:LJpot}
136   \end{equation}
137   and
138   \begin{equation}
139 < EQ Here
139 > V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
140 >        \boldsymbol{\Omega}_{j}) = \frac{|\mu_i||\mu_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[
141 >        \boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j}
142 >        -
143 >        3(\boldsymbol{\hat{u}}_i \cdot \mathbf{\hat{r}}_{ij}) %
144 >                (\boldsymbol{\hat{u}}_j \cdot \mathbf{\hat{r}}_{ij} )\biggr]
145   \label{lipidEq:dipolePot}
146   \end{equation}
147   Where $V_{\text{LJ}}$ is the Lennard-Jones potential and
# Line 44 | Line 150 | In $V_{\text{dipole}}$, $\mathbf{r}_{ij}$ is the vecto
150   parameters which scale the length and depth of the interaction
151   respectively, and $r_{ij}$ is the distance between beads $i$ and $j$.
152   In $V_{\text{dipole}}$, $\mathbf{r}_{ij}$ is the vector starting at
153 < bead$i$ and pointing towards bead $j$.  Vectors $\mathbf{\Omega}_i$
153 > bead $i$ and pointing towards bead $j$.  Vectors $\mathbf{\Omega}_i$
154   and $\mathbf{\Omega}_j$ are the orientational degrees of freedom for
155   beads $i$ and $j$.  $|\mu_i|$ is the magnitude of the dipole moment of
156   $i$, and $\boldsymbol{\hat{u}}_i$ is the standard unit orientation
157 < vector of $\boldsymbol{\Omega}_i$.
157 > vector rotated with Euler angles: $\boldsymbol{\Omega}_i$.
158  
159 < The model also allows for the bonded interactions of bonds, bends, and
160 < torsions.  The bonds between two beads on a chain are of fixed length,
161 < and are maintained according to the {\sc rattle} algorithm. \cite{fix}
162 < The bends are subject to a harmonic potential:
159 > The model also allows for the intra-molecular bend and torsion
160 > interactions.  The bond between two beads on a chain is of fixed
161 > length, and is maintained using the {\sc rattle}
162 > algorithm.\cite{andersen83} The bends are subject to a harmonic
163 > potential:
164   \begin{equation}
165 < eq here
165 > V_{\text{bend}}(\theta) = k_{\theta}( \theta - \theta_0 )^2
166   \label{lipidEq:bendPot}
167   \end{equation}
168 < where $fix$ scales the strength of the harmonic well, and $fix$ is the
169 < angle between bond vectors $fix$ and $fix$.  The torsion potential is
170 < given by:
168 > where $k_{\theta}$ scales the strength of the harmonic well, and
169 > $\theta$ is the angle between bond vectors
170 > (Fig.~\ref{lipidFig:lipidModel}). In addition, we have placed a
171 > ``ghost'' bend on the phospholipid head. The ghost bend is a bend
172 > potential which keeps the dipole roughly perpendicular to the
173 > molecular body, where $\theta$ is now the angle the dipole makes with
174 > respect to the {\sc head}-$\text{{\sc ch}}_2$ bond vector
175 > (Fig.~\ref{lipidFig:ghostBend}). This bend mimics the hinge between
176 > the phosphatidyl part of the PC head group and the remainder of the
177 > molecule.  The torsion potential is given by:
178   \begin{equation}
179 < eq here
179 > V_{\text{torsion}}(\phi) =  
180 >        k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0
181   \label{lipidEq:torsionPot}
182   \end{equation}
183   Here, the parameters $k_0$, $k_1$, $k_2$, and $k_3$ fit the cosine
184   power series to the desired torsion potential surface, and $\phi$ is
185 < the angle between bondvectors $fix$ and $fix$ along the vector $fix$
186 < (see Fig.:\ref{lipidFig:lipidModel}).  Long range interactions such as
187 < the Lennard-Jones potential are excluded for bead pairs involved in
188 < the same bond, bend, or torsion.  However, internal interactions not
189 < directly involved in a bonded pair are calculated.
185 > the angle the two end atoms have rotated about the middle bond
186 > (Fig.:\ref{lipidFig:lipidModel}).  Long range interactions such as the
187 > Lennard-Jones potential are excluded for atom pairs involved in the
188 > same bond, bend, or torsion.  However, long-range interactions for
189 > pairs of atoms not directly involved in a bond, bend, or torsion are
190 > calculated.
191  
192 + \begin{figure}
193 + \centering
194 + \includegraphics[width=0.5\linewidth]{ghostBendFig.eps}
195 + \caption[Depiction of the ``ghost'' bend]{The ``ghost'' bend is a bend potential added to restrain the motion of the dipole on the {\sc head} group. The potential follows Eq.~\ref{lipidEq:bendPot} where $\theta$ is now the angle that the dipole makes with the {\sc head}-$\text{{\sc ch}}_2$ bond vector.}
196 + \label{lipidFig:ghostBend}
197 + \end{figure}
198  
199 + All simulations presented here use a two-chain lipid as pictured in
200 + Fig.~\ref{lipidFig:lipidModel}.  The chains are both eight beads long,
201 + and their mass and Lennard Jones parameters are summarized in
202 + Table~\ref{lipidTable:tcLJParams}. The magnitude of the dipole moment
203 + for the head bead is 10.6~Debye (approximately half the magnitude of
204 + the dipole on the PC head group\cite{Cevc87}), and the bend and
205 + torsion parameters are summarized in
206 + Table~\ref{lipidTable:tcBendParams} and
207 + \ref{lipidTable:tcTorsionParams}.
208 +
209 + \begin{table}
210 + \caption[Lennard-Jones parameters for the two chain phospholipids]{THE LENNARD JONES PARAMETERS FOR THE TWO CHAIN PHOSPHOLIPIDS}
211 + \label{lipidTable:tcLJParams}
212 + \begin{center}
213 + \begin{tabular}{|l|c|c|c|c|}
214 + \hline
215 +     & mass (amu) & $\sigma$($\mbox{\AA}$)  & $\epsilon$ (kcal/mol) %
216 +        & $|\mathbf{\hat{\mu}}|$ (Debye) \\ \hline
217 + {\sc head} & 72  & 4.0 & 0.185 & 10.6 \\ \hline
218 + {\sc ch}\cite{Siepmann1998} & 13.02 & 4.0 & 0.0189 & 0.0 \\ \hline
219 + $\text{{\sc ch}}_2$\cite{Siepmann1998} & 14.03 & 3.95 & 0.18 & 0.0 \\ \hline
220 + $\text{{\sc ch}}_3$\cite{Siepmann1998} & 15.04 & 3.75 & 0.25 & 0.0 \\ \hline
221 + {\sc ssd}\cite{liu96:new_model} & 18.03 & 3.051 & 0.152 & 0.0 \\ \hline
222 + \end{tabular}
223 + \end{center}
224 + \end{table}
225 +
226 + \begin{table}
227 + \caption[Bend Parameters for the two chain phospholipids]{BEND PARAMETERS FOR THE TWO CHAIN PHOSPHOLIPIDS}
228 + \label{lipidTable:tcBendParams}
229 + \begin{center}
230 + \begin{tabular}{|l|c|c|}
231 + \hline
232 +   & $k_{\theta}$ ( kcal/($\text{mol deg}^2$) ) & $\theta_0$ ( deg ) \\ \hline
233 + {\sc ghost}-{\sc head}-$\text{{\sc ch}}_2$ & 0.00177 & 129.78 \\ \hline
234 + $x$-{\sc ch}-$y$ & 58.84 & 112.0 \\ \hline
235 + $x$-$\text{{\sc ch}}_2$-$y$ & 58.84 & 114.0 \\ \hline
236 + \end{tabular}
237 + \begin{minipage}{\linewidth}
238 + \begin{center}
239 + \vspace{2mm}
240 + All alkane parameters are based off of those from TraPPE.\cite{Siepmann1998}
241 + \end{center}
242 + \end{minipage}
243 + \end{center}
244 + \end{table}
245 +
246 + \begin{table}
247 + \caption[Torsion Parameters for the two chain phospholipids]{TORSION PARAMETERS FOR THE TWO CHAIN PHOSPHOLIPIDS}
248 + \label{lipidTable:tcTorsionParams}
249 + \begin{center}
250 + \begin{tabular}{|l|c|c|c|c|}
251 + \hline
252 + All are in kcal/mol $\rightarrow$ & $k_3$ & $k_2$ & $k_1$ & $k_0$ \\ \hline
253 + $x$-{\sc ch}-$y$-$z$ & 3.3254 & -0.4215 & -1.686 & 1.1661 \\ \hline
254 + $x$-$\text{{\sc ch}}_2$-$\text{{\sc ch}}_2$-$y$ & 5.9602 & -0.568 & -3.802 & 2.1586 \\ \hline
255 + \end{tabular}
256 + \begin{minipage}{\linewidth}
257 + \begin{center}
258 + \vspace{2mm}
259 + All alkane parameters are based off of those from TraPPE.\cite{Siepmann1998}
260 + \end{center}
261 + \end{minipage}
262 + \end{center}
263 + \end{table}
264 +
265 +
266 + \section{\label{lipidSec:furtherMethod}Further Methodology}
267 +
268 + As mentioned previously, the water model used throughout these
269 + simulations was the {\sc ssd/e} model of Fennell and
270 + Gezelter,\cite{fennell04} earlier forms of this model can be found in
271 + Ichiye \emph{et
272 + al}.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md} A
273 + discussion of the model can be found in Sec.~\ref{oopseSec:SSD}. As
274 + for the integration of the equations of motion, all simulations were
275 + performed in an orthorhombic periodic box with a thermostat on
276 + velocities, and an independent barostat on each Cartesian axis $x$,
277 + $y$, and $z$.  This is the $\text{NPT}_{xyz}$. integrator described in
278 + Sec.~\ref{oopseSec:integrate}. With $\tau_B = 1.5$~ps and $\tau_T =
279 + 1.2$~ps, the box volume stabilizes after 50~ps, and fluctuates about
280 + its equilibrium value by $\sim 0.6\%$, temperature fluctuations are
281 + about $\sim 1.4\%$ of their set value, and pressure fluctuations are
282 + the largest, varying as much as $\pm 250$~atm. However, such large
283 + fluctuations in pressure are typical for liquid state simulations.
284 +
285 +
286 + \subsection{\label{lipidSec:ExpSetup}Experimental Setup}
287 +
288 + Two main classes of starting configurations were used in this research:
289 + random and ordered bilayers.  The ordered bilayer simulations were all
290 + started from an equilibrated bilayer configuration at 300~K. The original
291 + configuration for the first 300~K run was assembled by placing the
292 + phospholipids centers of mass on a planar hexagonal lattice.  The
293 + lipids were oriented with their principal axis perpendicular to the plane.
294 + The bottom leaf simply mirrored the top leaf, and the appropriate
295 + number of water molecules were then added above and below the bilayer.
296 +
297 + The random configurations took more work to generate.  To begin, a
298 + test lipid was placed in a simulation box already containing water at
299 + the intended density.  The water molecules were then tested against
300 + the lipid using a 5.0~$\mbox{\AA}$ overlap test with any atom in the
301 + lipid.  This gave an estimate for the number of water molecules each
302 + lipid would displace in a simulation box. A target number of water
303 + molecules was then defined which included the number of water
304 + molecules each lipid would displace, the number of water molecules
305 + desired to solvate each lipid, and a factor to pad the initial box
306 + with a few extra water molecules.
307 +
308 + Next, a cubic simulation box was created that contained at least the
309 + target number of water molecules in an FCC lattice (the lattice was for ease of
310 + placement).  What followed was a RSA simulation similar to those of
311 + Chapt.~\ref{chapt:RSA}. The lipids were sequentially given a random
312 + position and orientation within the box.  If a lipid's position caused
313 + atomic overlap with any previously placed lipid, its position and
314 + orientation were rejected, and a new random placement site was
315 + attempted. The RSA simulation proceeded until all phospholipids had
316 + been placed.  After placement of all lipid molecules, water
317 + molecules with locations that overlapped with the atomic coordinates
318 + of the lipids were removed.
319 +
320 + Finally, water molecules were removed at random until the desired water
321 + to lipid ratio was achieved.  The typical low final density for these
322 + initial configurations was not a problem, as the box shrinks to an
323 + appropriate size within the first 50~ps of a simulation under the
324 + NPTxyz integrator.
325 +
326 + \subsection{\label{lipidSec:Configs}Configurations}
327 +
328 + The first class of simulations were started from ordered bilayers. All
329 + configurations consisted of 60 lipid molecules with 30 lipids on each
330 + leaf, and were hydrated with 1620 {\sc ssd/e} molecules. The original
331 + configuration was assembled according to Sec.~\ref{lipidSec:ExpSetup}
332 + and simulated for a length of 10~ns at 300~K. The other temperature
333 + runs were started from a configuration 7~ns in to the 300~K
334 + simulation. Their temperatures were modified with the thermostatting
335 + algorithm in the NPTxyz integrator. All of the temperature variants
336 + were also run for 10~ns, with only the last 5~ns being used for
337 + accumulation of statistics.
338 +
339 + The second class of simulations were two configurations started from
340 + randomly dispersed lipids in a ``gas'' of water. The first
341 + ($\text{R}_{\text{I}}$) was a simulation containing 72 lipids with
342 + 1800 {\sc ssd/e} molecules simulated at 300~K. The second
343 + ($\text{R}_{\text{II}}$) was 90 lipids with 1350 {\sc ssd/e} molecules
344 + simulated at 350~K. Both simulations were integrated for more than
345 + 20~ns to observe whether our model is capable of spontaneous
346 + aggregation into known phospholipid macro-structures:
347 + $\text{R}_{\text{I}}$ into a bilayer, and $\text{R}_{\text{II}}$ into
348 + a inverted rod.
349 +
350 + \section{\label{lipidSec:resultsDis}Results and Discussion}
351 +
352 + \subsection{\label{lipidSec:densProf}Density Profile}
353 +
354 + Fig.~\ref{lipidFig:densityProfile} illustrates the densities of the
355 + atoms in the bilayer systems normalized by the bulk density as a
356 + function of distance from the center of the box. The profile is taken
357 + along the bilayer normal (in this case the $z$ axis). The profile at
358 + 270~K shows several structural features that are largely smoothed out
359 + at 300~K. The left peak for the {\sc head} atoms is split at 270~K,
360 + implying that some freezing of the structure into a gel phase might
361 + already be occurring at this temperature. However, movies of the
362 + trajectories at this temperature show that the tails are very fluid,
363 + and have not gelled. But this profile could indicate that a phase
364 + transition may simply be beyond the time length of the current
365 + simulation, and that given more time the system may tend towards a gel
366 + phase. In all profiles, the water penetrates almost 5~$\mbox{\AA}$
367 + into the bilayer, completely solvating the {\sc head} atoms. The
368 + $\text{{\sc ch}}_3$ atoms, although mainly centered at the middle of
369 + the bilayer, show appreciable penetration into the head group
370 + region. This indicates that the chains have enough flexibility to bend
371 + back upward to allow the ends to explore areas around the {\sc head}
372 + atoms. It is unlikely that this is penetration from a lipid of the
373 + opposite face, as the lipids are only 12~$\mbox{\AA}$ in length, and
374 + the typical leaf spacing as measured from the {\sc head-head} spacing
375 + in the profile is 17.5~$\mbox{\AA}$.
376 +
377 + \begin{figure}
378 + \centering
379 + \includegraphics[width=\linewidth]{densityProfile.eps}
380 + \caption[The density profile of the lipid bilayers]{The density profile of the lipid bilayers along the bilayer normal. The black lines are the {\sc head} atoms, red lines are the {\sc ch} atoms, green lines are the $\text{{\sc ch}}_2$ atoms, blue lines are the $\text{{\sc ch}}_3$ atoms, and the magenta lines are the {\sc ssd} atoms.}
381 + \label{lipidFig:densityProfile}
382 + \end{figure}
383 +
384 +
385 + \subsection{\label{lipidSec:scd}$\text{S}_{\text{{\sc cd}}}$ Order Parameters}
386 +
387 + The $\text{S}_{\text{{\sc cd}}}$ order parameter is often reported in
388 + the experimental characterizations of phospholipids. It is obtained
389 + through deuterium NMR, and measures the ordering of the carbon
390 + deuterium bond in relation to the bilayer normal at various points
391 + along the chains. In our model, there are no explicit hydrogens, but
392 + the order parameter can be written in terms of the carbon ordering at
393 + each point in the chain:\cite{egberts88}
394 + \begin{equation}
395 + S_{\text{{\sc cd}}}  = \frac{2}{3}S_{xx} + \frac{1}{3}S_{yy}
396 + \label{lipidEq:scd1}
397 + \end{equation}
398 + Where $S_{ij}$ is given by:
399 + \begin{equation}
400 + S_{ij} = \frac{1}{2}\Bigl\langle(3\cos\Theta_i\cos\Theta_j
401 +        - \delta_{ij})\Bigr\rangle
402 + \label{lipidEq:scd2}
403 + \end{equation}
404 + Here, $\Theta_i$ is the angle the $i$th axis in the reference frame of
405 + the carbon atom makes with the bilayer normal. The brackets denote an
406 + average over time and molecules. The carbon atom axes are defined:
407 + \begin{itemize}
408 + \item $\mathbf{\hat{z}}\rightarrow$ vector from $C_{n-1}$ to $C_{n+1}$
409 + \item $\mathbf{\hat{y}}\rightarrow$ vector that is perpendicular to $z$ and
410 + in the plane through $C_{n-1}$, $C_{n}$, and $C_{n+1}$
411 + \item $\mathbf{\hat{x}}\rightarrow$ vector perpendicular to
412 + $\mathbf{\hat{y}}$ and $\mathbf{\hat{z}}$.
413 + \end{itemize}
414 + This assumes that the hydrogen atoms are always in a plane
415 + perpendicular to the $C_{n-1}-C_{n}-C_{n+1}$ plane.
416 +
417 + The order parameter has a range of $[1,-\frac{1}{2}]$. A value of 1
418 + implies full order aligned to the bilayer axis, 0 implies full
419 + disorder, and $-\frac{1}{2}$ implies full order perpendicular to the
420 + bilayer axis. The {\sc cd} bond vector for carbons near the head group
421 + are usually ordered perpendicular to the bilayer normal, with tails
422 + farther away tending toward disorder. This makes the order parameter
423 + negative for most carbons, and as such $|S_{\text{{\sc cd}}}|$ is more
424 + commonly reported than $S_{\text{{\sc cd}}}$.
425 +
426 + Fig.~\ref{lipidFig:scdFig} shows the $S_{\text{{\sc cd}}}$ order
427 + parameters for the bilayer system at 300~K. There is no appreciable
428 + difference in the plots for the various temperatures, however, there
429 + is a larger difference between our model's ordering, and the
430 + experimentally observed ordering of DMPC. As our values are closer to
431 + $-\frac{1}{2}$, this implies more ordering perpendicular to the normal
432 + than in a real system. This is due to the model having only one carbon
433 + group separating the chains from the top of the lipid. In DMPC, with
434 + the flexibility inherent in a multiple atom head group, as well as a
435 + glycerol linkage between the head group and the acyl chains, there is
436 + more loss of ordering by the point when the chains start.
437 +
438 + \begin{figure}
439 + \centering
440 + \includegraphics[width=\linewidth]{scdFig.eps}
441 + \caption[$\text{S}_{\text{{\sc cd}}}$ order parameter for our model]{A comparison of $|\text{S}_{\text{{\sc cd}}}|$ between our model (blue) and DMPC\cite{petrache00} (black) near 300~K.}
442 + \label{lipidFig:scdFig}
443 + \end{figure}
444 +
445 + \subsection{\label{lipidSec:p2Order}$P_2$ Order Parameter}
446 +
447 + The $P_2$ order parameter allows us to measure the amount of
448 + directional ordering that exists in the bodies of the molecules making
449 + up the bilayer. Each lipid molecule can be thought of as a cylindrical
450 + rod with the head group at the top. If all of the rods are perfectly
451 + aligned, the $P_2$ order parameter will be $1.0$. If the rods are
452 + completely disordered, the $P_2$ order parameter will be 0. For a
453 + collection of unit vectors pointing along the principal axes of the
454 + rods, the $P_2$ order parameter can be solved via the following
455 + method.\cite{zannoni94}
456 +
457 + Define an ordering tensor $\overleftrightarrow{\mathsf{Q}}$, such that,
458 + \begin{equation}
459 + \overleftrightarrow{\mathsf{Q}} = \frac{1}{N}\sum_i^N %
460 +        \begin{pmatrix} %
461 +        u_{ix}u_{ix}-\frac{1}{3} & u_{ix}u_{iy} & u_{ix}u_{iz} \\
462 +        u_{iy}u_{ix} & u_{iy}u_{iy}-\frac{1}{3} & u_{iy}u_{iz} \\
463 +        u_{iz}u_{ix} & u_{iz}u_{iy} & u_{iz}u_{iz}-\frac{1}{3} %
464 +        \end{pmatrix}
465 + \label{lipidEq:po1}
466 + \end{equation}
467 + Where the $u_{i\alpha}$ is the $\alpha$ element of the unit vector
468 + $\mathbf{\hat{u}}_i$, and the sum over $i$ averages over the whole
469 + collection of unit vectors. This allows the tensor to be written:
470 + \begin{equation}
471 + \overleftrightarrow{\mathsf{Q}} = \frac{1}{N}\sum_i^N \biggl[
472 +        \mathbf{\hat{u}}_i \otimes \mathbf{\hat{u}}_i
473 +        - \frac{1}{3} \cdot \mathsf{1} \biggr]
474 + \label{lipidEq:po2}
475 + \end{equation}
476 +
477 + After constructing the tensor, diagonalizing
478 + $\overleftrightarrow{\mathsf{Q}}$ yields three eigenvalues and
479 + eigenvectors. The eigenvector associated with the largest eigenvalue,
480 + $\lambda_{\text{max}}$, is the director axis  for the system of unit
481 + vectors. The director axis is the average direction all of the unit vectors
482 + are pointing. The $P_2$ order parameter is then simply
483 + \begin{equation}
484 + \langle P_2 \rangle = \frac{3}{2}\lambda_{\text{max}}
485 + \label{lipidEq:po3}
486 + \end{equation}
487 +
488 + Table~\ref{lipidTab:blSummary} summarizes the $P_2$ values for the
489 + bilayers, as well as the dipole orientations. The unit vector for the
490 + lipid molecules was defined by finding the moment of inertia for each
491 + lipid, then setting $\mathbf{\hat{u}}$ to point along the axis of
492 + minimum inertia. For the {\sc head} atoms, the unit vector simply
493 + pointed in the same direction as the dipole moment. For the lipid
494 + molecules, the ordering was consistent across all temperatures, with
495 + the director pointed along the $z$ axis of the box. More
496 + interestingly, is the high degree of ordering the dipoles impose on
497 + the {\sc head} atoms. The directors for the dipoles themselves
498 + consistently pointed along the plane of the bilayer, with the
499 + directors anti-aligned on the top and bottom leaf.
500 +
501 + \begin{table}
502 + \caption[Structural properties of the bilayers]{BILAYER STRUCTURAL PROPERTIES AS A FUNCTION OF TEMPERATURE}
503 + \label{lipidTab:blSummary}
504 + \begin{center}
505 + \begin{tabular}{|c|c|c|c|c|}
506 + \hline
507 + Temperature (K) & $\langle L_{\perp}\rangle$ ($\mbox{\AA}$) & %
508 +        $\langle A_{\parallel}\rangle$ ($\mbox{\AA}^2$) & %
509 +        $\langle P_2\rangle_{\text{Lipid}}$ & %
510 +        $\langle P_2\rangle_{\text{{\sc head}}}$ \\ \hline
511 + 270 & 18.1 & 58.1 & 0.253 & 0.494 \\ \hline
512 + 275 & 17.2 & 56.7 & 0.295 & 0.514 \\ \hline
513 + 277 & 16.9 & 58.0 & 0.301 & 0.541 \\ \hline
514 + 280 & 17.4 & 58.0 & 0.274 & 0.488 \\ \hline
515 + 285 & 16.9 & 57.6 & 0.270 & 0.616 \\ \hline
516 + 290 & 17.0 & 57.6 & 0.263 & 0.534 \\ \hline
517 + 293 & 17.5 & 58.0 & 0.227 & 0.643 \\ \hline
518 + 300 & 16.9 & 57.6 & 0.315 & 0.536 \\ \hline
519 + \end{tabular}
520 + \end{center}
521 + \end{table}
522 +
523 + \subsection{\label{lipidSec:miscData}Further Structural Data}
524 +
525 + Also summarized in Table~\ref{lipidTab:blSummary}, are the bilayer
526 + thickness ($\langle L_{\perp}\rangle$) and area per lipid ($\langle
527 + A_{\parallel}\rangle$). The bilayer thickness was measured from the
528 + peak to peak {\sc head} atom distance in the density profiles. The
529 + area per lipid data compares favorably with values typically seen for
530 + DMPC (60.0~$\mbox{\AA}^2$ at 303~K)\cite{petrache00}. Although our
531 + values are lower this is most likely due to the shorter chain length
532 + of our model (8 versus 14 for DMPC).
533 +
534 + \subsection{\label{lipidSec:diffusion}Lateral Diffusion Constants}
535 +
536 + The lateral diffusion constant, $D_L$, is the constant characterizing
537 + the diffusive motion of the lipid molecules within the plane of the bilayer. It
538 + is given by the following Einstein relation:\cite{allen87:csl}
539 + \begin{equation}
540 + D_L = \lim_{t\rightarrow\infty}\frac{1}{4t}\langle |\mathbf{r}(t)
541 +        - \mathbf{r}(0)|^2\rangle
542 + \end{equation}
543 + Where $\mathbf{r}(t)$ is the $xy$ position of the lipid at time $t$
544 + (assuming the $z$-axis is parallel to the bilayer normal).
545 +
546 + Fig.~\ref{lipidFig:diffusionFig} shows the lateral diffusion constants
547 + as a function of temperature. There is a definite increase in the
548 + lateral diffusion with higher temperatures, which is exactly what one
549 + would expect with greater fluidity of the chains. However, the
550 + diffusion constants are two orders of magnitude smaller than those
551 + typical of DPPC.\cite{Cevc87} This is counter-intuitive as the DPPC
552 + molecule is sterically larger and heavier than our model. This could
553 + be an indication that our model's chains are too interwoven and hinder
554 + the motion of the lipid or that the dipolar head groups are too
555 + tightly bound to each other. In contrast, the diffusion constant of
556 + the {\sc ssd} water, $9.84\times 10^{-6}\,\text{cm}^2/\text{s}$, is
557 + reasonably close to the bulk water diffusion constant ($2.2999\times
558 + 10^{-5}\,\text{cm}^2/\text{s}$).\cite{Holz00}
559 +
560 + \begin{figure}
561 + \centering
562 + \includegraphics[width=\linewidth]{diffusionFig.eps}
563 + \caption[The lateral diffusion constants versus temperature]{The lateral diffusion constants for the bilayers as a function of temperature.}
564 + \label{lipidFig:diffusionFig}
565 + \end{figure}
566 +
567 + \subsection{\label{lipidSec:randBilayer}Bilayer Aggregation}
568 +
569 + A very important accomplishment for our model is its ability to
570 + spontaneously form bilayers from a randomly dispersed starting
571 + configuration. Fig.~\ref{lipidFig:blImage} shows an image sequence for
572 + the bilayer aggregation. After 3.0~ns, the basic form of the bilayer
573 + can already be seen. By 7.0~ns, the bilayer has a lipid bridge
574 + stretched across the simulation box to itself that will turn out to be
575 + very long lived ($\sim$20~ns), as well as a water pore, that will
576 + persist for the length of the current simulation. At 24~ns, the lipid
577 + bridge has broken, and the bilayer is still integrating the lipid
578 + molecules from the bridge into itself. However, the water pore is
579 + still present at 24~ns.
580 +
581 + \begin{figure}
582 + \centering
583 + \includegraphics[width=\linewidth]{bLayerImage.eps}
584 + \caption[Image sequence of the bilayer aggregation]{Image sequence of the bilayer aggregation. The blue beads are the {\sc head} atoms the grey beads are the chains, and the red and white bead are the water molecules. A box has been drawn around the periodic image.}
585 + \label{lipidFig:blImage}
586 + \end{figure}
587 +
588 + \subsection{\label{lipidSec:randIrod}Inverted Rod Aggregation}
589 +
590 + Fig.~\ref{lipidFig:iRimage} shows a second aggregation sequence
591 + simulated in this research. Here the fraction of water had been
592 + significantly decreased to observe how the model would respond. After
593 + 1.5~ns, The main body of water in the system has already collected
594 + into a central water channel. By 10.0~ns, the channel has widened
595 + slightly, but there are still many water molecules permeating the
596 + lipid macro-structure. At 23.0~ns, the central water channel has
597 + stabilized and several smaller water channels have been absorbed by
598 + the main one. However, there is still an appreciable water
599 + concentration throughout the lipid structure.
600 +
601 + \begin{figure}
602 + \centering
603 + \includegraphics[width=\linewidth]{iRodImage.eps}
604 + \caption[Image sequence of the inverted rod aggregation]{Image sequence of the inverted rod aggregation. color scheme is the same as in Fig.~\ref{lipidFig:blImage}.}
605 + \label{lipidFig:iRimage}
606 + \end{figure}
607 +
608 + \section{\label{lipidSec:Conclusion}Conclusion}
609 +
610 + We have presented a simple unified-atom phospholipid model capable of
611 + spontaneous aggregation into a bilayer and an inverted rod
612 + structure. The time scales of the macro-molecular aggregations are
613 + approximately 24~ns. In addition the model's properties have been
614 + explored over a range of temperatures through prefabricated
615 + bilayers. No freezing transition is seen in the temperature range of
616 + our current simulations. However, structural information from 270~K
617 + may imply that a freezing event is on a much longer time scale than
618 + that explored in this current research. Further studies of this system
619 + could extend the time length of the simulations at the low
620 + temperatures to observe whether lipid crystallization can occur within
621 + the framework of this model.
622 +
623 + Potential problems that may be obstacles in further research, is the
624 + lack of detail in the head region. As the chains are almost directly
625 + attached to the {\sc head} atom, there is no buffer between the
626 + actions of the head group and the tails. Another disadvantage of the
627 + model is the dipole approximation will alter results when details
628 + concerning a charged solute's interactions with the bilayer. However,
629 + it is important to keep in mind that the dipole approximation can be
630 + kept an advantage by examining solutes that do not require point
631 + charges, or at the least, require only dipole approximations
632 + themselves. Other advantages of the model include the ability to alter
633 + the size of the unified-atoms so that the size of the lipid can be
634 + increased without adding to the number of interactions in the
635 + system. However, what sets our model apart from other current
636 + simplified models,\cite{goetz98,marrink04} is the information gained
637 + by observing the ordering of the head groups dipole's in relation to
638 + each other and the solvent without the need for point charges and the
639 + Ewald sum.

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines