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 1087 by mmeineke, Fri Mar 5 22:16:34 2004 UTC vs.
Revision 1089 by mmeineke, Mon Mar 8 22:15:54 2004 UTC

# Line 7 | Line 7 | 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 phase ($L_{\beta}$) of
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
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
# Line 21 | Line 21 | computationally expensive Ewald sum or the faster, par
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 Ewald sum or the faster, particle mesh Ewald
25 < sum.\cite{nina:2002,norberg:2000,patra:2003} These factors all limit
26 < the system size and time scales of bilayer simulations.
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 well beyond the range of current simulation technology. One
30 < such example is the observance of a ripple phase
31 < ($P_{\beta^{\prime}}$) between the $L_{\beta}$ and $L_{\alpha}$ phases
32 < of certain phospholipid bilayers.\cite{katsaras00,sengupta00} These
33 < ripples are shown to have periodicity on the order of
34 < 100-200~$\mbox{\AA}$. A simulation on this length scale would have
35 < approximately 1,300 lipid molecules with an additional 25 water
36 < molecules per lipid to fully solvate the bilayer. A simulation of this
37 < size is impractical with current atomistic models.
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   The time and length scale limitations are most striking in transport
49 < phenomena.  Due to the fluid-like properties of a lipid membrane, not
50 < all diffusion across the membrane happens at pores.  Some molecules of
51 < interest may incorporate themselves directly into the membrane.  Once
52 < here, they may possess an appreciable waiting time (on the order of
53 < 10's to 100's of nanoseconds) within the bilayer. Such long simulation
54 < times are nearly impossible to obtain when integrating the system with
55 < atomistic detail.
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 Liposky\cite{goetz98} is to model the entire
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}$), and eliminates
81 < the 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
90 < molecules.
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
# Line 114 | Line 124 | The long range interactions between lipids are given b
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 130 | Line 141 | V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_
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]
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 145 | Line 156 | The model also allows for the bonded interactions bend
156   $i$, and $\boldsymbol{\hat{u}}_i$ is the standard unit orientation
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 156 | 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
175 < (Fig.~\ref{lipidFig:ghostBend}). The torsion potential is given 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 171 | Line 185 | same bond, bend, or torsion.  However, internal intera
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   \begin{figure}
193   \centering
194 < \includegraphics[width=\linewidth]{ghostBendFig.eps}
195 < \caption[Depiction of the ``ghost'' bend]{The ``ghost'' bend is a bend potential added to constrain 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.}
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 chained lipid as pictured in
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{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 216 | 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 229 | 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:integrate}.
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
259 < 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 factor to pad the
305 < initial box with a few extra water molecules.
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 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 placement of all lipid molecules, water
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
321 < number of waters per lipid was reached.  The typical low final density
322 < for these initial configurations was not a problem, as the box shrinks
323 < to an appropriate size within the first 50~ps of a simulation in the
324 < $\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}Configurations}
327  
328 < The first class of simulations were started from ordered
329 < bilayers. They were all configurations consisting of 60 lipid
330 < molecules with 30 lipids on each leaf, and were hydrated with 1620
331 < {\sc ssd} molecules. The original configuration was assembled
332 < according to Sec.~\ref{lipidSec:ExpSetup} and simulated for a length
333 < of 10~ns at 300~K. The other temperature runs were started from a
334 < frame 7~ns into the 300~K simulation. Their temperatures were reset
335 < with the thermostating algorithm in the $\text{NPT}_{xyz}$
336 < integrator. All of the temperature variants were also run for 10~ns,
337 < with only the last 5~ns being used for accumulation of statistics.
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} molecules simulated at 300~K. The second
343 < ($\text{R}_{\text{II}}$) was 90 lipids with 1350 {\sc ssd} molecules
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, and illustrate the spontaneous aggregation of the lipid model
346 < into phospholipid macro-structures: $\text{R}_{\text{I}}$ into a
347 < bilayer, and $\text{R}_{\text{II}}$ into a inverted rod.
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  
314 \subsection{\label{lipidSec:diffusion}Lateral Diffusion Constants}
315
316 The lateral diffusion constant, $D_L$, is the constant characterizing
317 the diffusive motion of the lipid within the plane of the bilayer. It
318 is given by the following Einstein relation valid at long
319 times:\cite{allen87:csl}
320 \begin{equation}
321 2tD_L = \frac{1}{2}\langle |\mathbf{r}(t) - \mathbf{r}(0)|^2\rangle
322 \end{equation}
323 Where $\mathbf{r}(t)$ is the position of the lipid at time $t$, and is
324 constrained to lie within a plane. For the bilayer simulations the
325 plane of constrained motion was that perpendicular to the bilayer
326 normal, namely the $xy$-plane.
327
328 Fig.~\ref{lipidFig:diffusionFig} shows the lateral diffusion constants
329 as a function of temperature. There is a definite increase in the
330 lateral diffusion with higher temperatures, which is exactly what one
331 would expect with greater fluidity of the chains. However, the
332 diffusion constants are all two orders of magnitude smaller than those
333 typical of DPPC.\cite{Cevc87} This is counter-intuitive as the DPPC
334 molecule is sterically larger and heavier than our model. This could
335 be an indication that our model's chains are too interwoven and hinder
336 the motion of the lipid, or that a simplification in the model's
337 forces has led to a slowing of diffusive behavior within the
338 bilayer. In contrast, the diffusion constant of the {\sc ssd} water,
339 $9.84\times 10^{-6}\,\text{cm}^2/\text{s}$, compares favorably with
340 that of bulk water ($2.2999\times
341 10^{-5}\,\text{cm}^2/\text{s}$\cite{Holz00}).
342
343 \begin{figure}
344 \centering
345 \includegraphics[width=\linewidth]{diffusionFig.eps}
346 \caption[The lateral diffusion constants versus temperature]{The lateral diffusion constants for the bilayers as a function of temperature.}
347 \label{lipidFig:diffusionFig}
348 \end{figure}
349
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
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 < by 300~K. The left peak for the {\sc head} atoms is split at 270~K,
360 < implying that some freezing of the structure might already be occurring
361 < at this temperature. From the dynamics, the tails at this temperature
362 < are very much fluid, but the profile could indicate that a phase
363 < transition may simply be beyond the length scale of the current
364 < simulation. In all profiles, the water penetrates almost
365 < 5~$\mbox{\AA}$ into the bilayer, completely solvating the {\sc head}
366 < atoms. The $\text{{\sc ch}}_3$ atoms although mainly centered at the
367 < middle of the bilayer, show appreciable penetration into the head
368 < group region. This indicates that the chains have enough mobility to
369 < bend back upward to allow the ends to explore areas around the {\sc
370 < head} atoms. It is unlikely that this is penetration from a lipid of
371 < the opposite face, as the lipids are only 12~$\mbox{\AA}$ in length,
372 < and the typical leaf spacing as measured from the {\sc head-head}
373 < spacing in the profile is 17.5~$\mbox{\AA}$.
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
# Line 400 | Line 404 | $\mathbf{\hat{z}}\rightarrow$ vector from $C_{n-1}$ to
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 < $\mathbf{\hat{z}}\rightarrow$ vector from $C_{n-1}$ to $C_{n+1}$;
408 < $\mathbf{\hat{y}}\rightarrow$ vector that is perpendicular to $z$ and
409 < in the plane through $C_{n-1}$, $C_{n}$, and $C_{n+1}$;
410 < $\mathbf{\hat{x}}\rightarrow$ vector perpendicular to
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
# Line 418 | Line 426 | is a larger difference between our models ordering, an
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 models ordering, and that of
430 < DMPC. As our values are closer to $-\frac{1}{2}$, this implies more
431 < ordering perpendicular to the normal than in a real system. This is
432 < due to the model having only one carbon group separating the chains
433 < from the top of the lipid. In DMPC, with the flexibility inherent in a
434 < multiple atom head group, as well as a glycerol linkage between the
435 < head group and the acyl chains, there is more loss of ordering by the
436 < point when the chains start.
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
# Line 437 | Line 445 | directional ordering that exists in the bilayer. Each
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 bilayer. Each lipid molecule
449 < can be thought of as a cylindrical tube with the head group at the
450 < top. If all of the cylinders are perfectly aligned, the $P_2$ order
451 < parameter will be $1.0$. If the cylinders are completely dispersed,
452 < the $P_2$ order parameter will be 0. For a collection of unit vectors,
453 < the $P_2$ order parameter can be solved via the following
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 matrix $\mathbf{Q}$, such that,
457 > Define an ordering tensor $\overleftrightarrow{\mathsf{Q}}$, such that,
458   \begin{equation}
459 < \mathbf{Q} = \frac{1}{N}\sum_i^N %
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} \\
# Line 457 | Line 466 | collection of unit vectors. This allows the matrix ele
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 matrix element
461 < $Q_{\alpha\beta}$ to be written:
469 > collection of unit vectors. This allows the tensor to be written:
470   \begin{equation}
471 < Q_{\alpha\beta} = \langle u_{\alpha}u_{\beta} -
472 <        \frac{1}{3}\delta_{\alpha\beta} \rangle
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 < Having constructed the matrix, diagonalizing $\mathbf{Q}$ yields three
478 < eigenvalues and eigenvectors. The eigenvector associated with the
479 < largest eigenvalue, $\lambda_{\text{max}}$, is the director for the
480 < system of unit vectors. The director is the average direction all of
481 < the unit vectors are pointing. The $P_2$ order parameter is then
482 < simply
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}
# Line 485 | Line 494 | the {\sc head} atoms. The directors for the dipoles co
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 consistently
498 < pointed along the plane of the bilayer, with the directors
499 < anti-aligned on the top and bottom leaf.
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.}
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|}
# Line 511 | Line 520 | Temperature (K) & $\langle L_{\perp}\rangle$ ($\mbox{\
520   \end{center}
521   \end{table}
522  
523 < \subsection{\label{lipidSec:miscData}Further Bilayer Data}
523 > \subsection{\label{lipidSec:miscData}Further Structural Data}
524  
525   Also summarized in Table~\ref{lipidTab:blSummary}, are the bilayer
526 < thickness and area per lipid. The bilayer thickness was measured from
527 < the peak to peak {\sc head} atom distance in the density profiles. The
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 are
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
# Line 531 | Line 574 | bridge is dispersed, and the bilayer is still integrat
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 is dispersed, and the bilayer is still integrating the lipid
578 < molecules from the bridge into itself, and has still been unable to
579 < disperse the water pore.
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
# Line 549 | Line 592 | slightly, but there are still many sub channels permea
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 sub channels permeating the lipid
596 < macro-structure. At 23.0~ns, the central water channel has stabilized
597 < and several smaller water channels have been absorbed to main
598 < one. However, there are still several other channels that persist
599 < through the lipid structure.
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
# Line 564 | Line 607 | We have presented a phospholipid model capable of spon
607  
608   \section{\label{lipidSec:Conclusion}Conclusion}
609  
610 < We have presented a phospholipid model capable of spontaneous
611 < aggregation into a bilayer and an inverted rod structure. The time
612 < scales of the macro-molecular aggregations are in excess of 24~ns. In
613 < addition the model's bilayer properties have been explored over a
614 < range of temperatures through prefabricated bilayers. No freezing
615 < transition is seen in the temperature range of our current
616 < simulations. However, structural information from the lowest
617 < temperature may imply that a freezing event is on a much longer time
618 < scale than that explored in this current research. Further studies of
619 < this system could extend the time length of the simulations at the low
620 < temperatures to observe whether lipid crystallization occurs within the
621 < framework of this model.
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