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 1061 by mmeineke, Fri Feb 20 17:52:57 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 < In the past 10 years, computer speeds have allowed for the atomistic
8 < simulation of phospholipid bilayers.  These simulations have ranged
9 < from simulation of the gel phase ($L_{\beta}$) of
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 64 to 256
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 < This is due to the expense of the computer calculations involved when
19 < performing these simulations.  To properly hydrate a bilayer, one
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 the computationally
24 < expensive Ewald sum or its slightly faster derivative particle mesh
25 < Ewald sum.\cite{nina:2002,norberg:2000,patra:2003} These factors all
26 < limit the potential size and time lengths of bilayer simulations.
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   Unfortunately, much of biological interest happens on time and length
29 < scales infeasible with current simulation. One such example is the
30 < observance of a ripple phase ($P_{\beta^{\prime}}$) between the
31 < $L_{\beta}$ and $L_{\alpha}$ phases of certain phospholipid
32 < bilayers.\cite{katsaras00,sengupta00} These ripples are shown to
33 < have periodicity on the order of 100-200~$\mbox{\AA}$. A simulation on
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 < Another class of simulations to consider, are those dealing with the
42 < diffusion of molecules through a bilayer.  Due to the fluid-like
43 < properties of a lipid membrane, not all diffusion across the membrane
44 < happens at pores.  Some molecules of interest may incorporate
45 < themselves directly into the membrane.  Once here, they may possess an
46 < appreciable waiting time (on the order of 10's to 100's of
44 < nanoseconds) within the bilayer. Such long simulation times are
45 < difficulty to obtain when integrating the system with atomistic
46 < detail.
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 < Addressing these issues, several schemes have been proposed.  One
49 < approach by Goetz and Liposky\cite{goetz98} is to model the entire
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 > 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 similar
63 < model proposed by Marrinck \emph{et. al}.\cite{marrink04}~uses a
64 < similar technique for modeling polar and non-polar interactions with
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. While the solvent spheres are kept charge-neutral and
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 simplifications over atomistic detail.  The model uses
74 < Lennard-Jones spheres for the head and tail groups of 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}$), eliminating the
81 < need for a costly Ewald sum.  
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
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 systems in this research mimic the entropic
88 < contribution to the hydrophobic effect due to hydrogen-bond network
89 < deformation around a non-polar entity, \emph{i.e.}~the phospholipid.
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
95 < used in these simulations.  As well as clarification about the water
96 < model and integration techniques.  The various simulation setups
97 < explored in this research are outlined in
98 < Sec.~\ref{lipidSec:ExpSetup}. Sec.~\ref{lipidSec:Results} and
99 < Sec.~\ref{lipidSec:Discussion} give a summary of the results and
100 < interpretation of those results respectively.  Finally, the
90 < conclusions of this chapter are presented in
91 < Sec.~\ref{lipidSec:Conclusion}.
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  
# Line 109 | Line 118 | given phospholipid head group.  The atoms of the tail
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 $\sigma_{\text{tail}}$ and $\epsilon_{\text{tail}}$.
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   V_{\text{LJ}}(r_{ij}) =
131          4\epsilon_{ij} \biggl[
# Line 129 | Line 140 | V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_
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 <        \frac{3(\boldsymbol{\hat{u}}_i \cdot \mathbf{r}_{ij}) %
144 <                (\boldsymbol{\hat{u}}_j \cdot \mathbf{r}_{ij}) }
134 <                {r^{2}_{ij}} \biggr]
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 140 | 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 bends, and torsions.
160 < The bond between two beads on a chain is of fixed length, and is
161 < maintained according to the {\sc rattle} algorithm.\cite{andersen83}
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   V_{\text{bend}}(\theta) = k_{\theta}( \theta - \theta_0 )^2
166   \label{lipidEq:bendPot}
# Line 157 | Line 168 | $\theta$ is the angle between bond vectors
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 adds a
172 < potential to keep the dipole pointed along the bilayer surface, where
173 < $theta$ is now the angle the dipole makes with respect to the {\sc
174 < head}-$\text{{\sc ch}}_2$ bond vector. The torsion potential is given
175 < by:
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   V_{\text{torsion}}(\phi) =  
180          k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0
# Line 172 | Line 185 | Lennard-Jones potential are excluded for atom pairs in
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, internal interactions not
189 < directly involved in a bonded pair are calculated.
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 < All simulations presented here use a two chained lipid as pictured in
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, and the bend and torsion parameters
204 < are summarized in Table~\ref{lipidTable:tcBendParams} and
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{The Lennard Jones Parameters for the two chain phospholipids.}
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|}
213 > \begin{tabular}{|l|c|c|c|c|}
214   \hline
215 <     & mass (amu) & $\sigma$($\mbox{\AA}$)  & $\epsilon$ (kcal/mol) \\ \hline
216 < {\sc head} & 72  & 4.0 & 0.185 \\ \hline
217 < {\sc ch}\cite{Siepmann1998} & 13.02 & 4.0 & 0.0189 \\ \hline
218 < $\text{{\sc ch}}_2$\cite{Siepmann1998} & 14.03 & 3.95 & 0.18 \\ \hline
219 < $\text{{\sc ch}}_3$\cite{Siepmann1998} & 15.04 & 3.75 & 0.25 \\ \hline
220 < {\sc ssd}\cite{liu96:new_model} & 18.03 & 3.051 & 0.152 \\ \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. All alkane parameters are based off of those from TraPPE.\cite{Siepmann1998}}
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|}
# Line 210 | Line 234 | $x$-$\text{{\sc ch}}_2$-$y$ & 58.84 & 114.0 \\ \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. Alkane parameters based on TraPPE.\cite{Siepmann1998}}
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|}
# Line 223 | Line 253 | $x$-$\text{{\sc ch}}_2$-$\text{{\sc ch}}_2$-$y$ & 5.96
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} model of
270 < Ichiye.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md} A
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}$. ensemble described in
278 < Sec.~\ref{oopseSec:Ensembles}.
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 starting configuration classes were used in this research:
289 < random and ordered bilayers.  The ordered bilayer starting
290 < configurations were all started from an equilibrated bilayer at
291 < 300~K. The original configuration for the first 300~K run was
292 < assembled by placing the phospholipids centers of mass on a planar
293 < hexagonal lattice.  The lipids were oriented with their long axis
294 < perpendicular to the plane.  The second leaf simply mirrored the first
295 < leaf, and the appropriate number of waters were then added above and
253 < below the bilayer.
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 waters were then tested for overlap with
300 < the lipid using a 5.0~$\mbox{\AA}$ buffer distance.  This gave an
301 < estimate for the number of waters each lipid would displace in a
302 < simulation box. A target number of waters was then defined which
303 < included the number of waters each lipid would displace, the number of
304 < waters desired to solvate each lipid, and a fudge factor to pad the
305 < initialization.
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 waters in an FCC lattice (the lattice was for ease of
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 adsorbed lipid, its position and
314 < orientation were rejected, and a new random adsorption site was
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 adsorbed.  After adsorption, all water molecules with locations
317 < that overlapped with the atomic coordinates of the lipids were
318 < removed.
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 one by one at random until the
321 < desired number of waters per lipid was reached.  The typical low final
322 < density for these initial configurations was not a problem, as the box
323 < would shrink to an appropriate size within the first 50~ps of a
324 < simulation in the $\text{NPT}_{xyz}$ ensemble.
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}The simulation configurations}
326 > \subsection{\label{lipidSec:Configs}Configurations}
327  
328 < Table ~\ref{lipidTable:simNames} summarizes the names and important
329 < details of the simulations.  The B set of simulations were all started
330 < in an ordered bilayer and observed over a period of 10~ns. Simulation
331 < RL was integrated for approximately 20~ns starting from a random
332 < configuration as an example of spontaneous bilayer aggregation.
333 < Lastly, simulation RH was also started from a random configuration,
334 < but with a lesser water content and higher temperature to show the
335 < spontaneous aggregation of an inverted hexagonal lamellar phase.
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