| 28 |
|
hydrophobic interior of the membrane, and for the hydrophobic tails |
| 29 |
|
to be exposed to the aqueous environment\cite{Sasaki2004}. A number |
| 30 |
|
of studies indicate that the flipping of phospholipids occurs |
| 31 |
< |
rapidly in the eukaryotic ER and the bacterial cytoplasmic membrane |
| 32 |
< |
via a bi-directional, facilitated diffusion process requiring no |
| 33 |
< |
metabolic energy input. Another system of interest would be the |
| 34 |
< |
distribution of sites occupied by inhaled anesthetics in membrane. |
| 35 |
< |
Although the physiological effects of anesthetics have been |
| 36 |
< |
extensively studied, the controversy over their effects on lipid |
| 37 |
< |
bilayers still continues. Recent deuterium NMR measurements on |
| 38 |
< |
halothane in POPC lipid bilayers suggest the anesthetics are |
| 39 |
< |
primarily located at the hydrocarbon chain region\cite{Baber1995}. |
| 40 |
< |
Infrared spectroscopy experiments suggest that halothane in DMPC |
| 41 |
< |
lipid bilayers lives near the membrane/water |
| 42 |
< |
interface\cite{Lieb1982}. |
| 31 |
> |
rapidly in the eukaryotic endoplasmic reticulum and the bacterial |
| 32 |
> |
cytoplasmic membrane via a bi-directional, facilitated diffusion |
| 33 |
> |
process requiring no metabolic energy input. Another system of |
| 34 |
> |
interest is the distribution of sites occupied by inhaled |
| 35 |
> |
anesthetics in membrane. Although the physiological effects of |
| 36 |
> |
anesthetics have been extensively studied, the controversy over |
| 37 |
> |
their effects on lipid bilayers still continues. Recent deuterium |
| 38 |
> |
NMR measurements on halothane on POPC lipid bilayers suggest the |
| 39 |
> |
anesthetics are primarily located at the hydrocarbon chain |
| 40 |
> |
region\cite{Baber1995}. However, infrared spectroscopy experiments |
| 41 |
> |
suggest that halothane in DMPC lipid bilayers lives near the |
| 42 |
> |
membrane/water interface\cite{Lieb1982}. |
| 43 |
|
|
| 44 |
|
Molecular dynamics simulations have proven to be a powerful tool for |
| 45 |
|
studying the functions of biological systems, providing structural, |
| 69 |
|
where $r_{ij}$ is the vector between molecule $i$ and molecule $j$, |
| 70 |
|
$\Omega _i$ and $\Omega _j$ are the orientational degrees of freedom |
| 71 |
|
for molecule $i$ and molecule $j$ respectively. |
| 72 |
< |
\[ |
| 73 |
< |
V_{LJ} (r_{ij} ) = 4\varepsilon _{ij} \left[ {\left( {\frac{{\sigma |
| 74 |
< |
_{ij} }}{{r_{ij} }}} \right)^{12} - \left( {\frac{{\sigma _{ij} |
| 75 |
< |
}}{{r_{ij} }}} \right)^6 } \right] |
| 76 |
< |
\] |
| 77 |
< |
\[ |
| 78 |
< |
V_{dp} (r_{ij} ,\Omega _i ,\Omega _j ) = \frac{1}{{4\pi \varepsilon |
| 79 |
< |
_0 }}\left[ {\frac{{\mu _i \cdot \mu _j }}{{r_{ij}^3 }} - |
| 80 |
< |
\frac{{3\left( {\mu _i \cdot r_{ij} } \right)\left( {\mu _i \cdot |
| 81 |
< |
r_{ij} } \right)}}{{r_{ij}^5 }}} \right] |
| 82 |
< |
\] |
| 83 |
< |
\[ |
| 84 |
< |
V_{sticky} (r_{ij} ,\Omega _i ,\Omega _j ) = v_0 [s(r_{ij} )w(r_{ij} |
| 85 |
< |
,\Omega _i ,\Omega _j ) + s'(r_{ij} )w'(r_{ij} ,\Omega _i ,\Omega _j |
| 86 |
< |
)] |
| 87 |
< |
\] |
| 72 |
> |
\begin{eqnarray*} |
| 73 |
> |
V_{LJ} (r_{ij} ) &= &4\varepsilon _{ij} \left[ {\left( |
| 74 |
> |
{\frac{{\sigma _{ij} }}{{r_{ij} }}} \right)^{12} - \left( |
| 75 |
> |
{\frac{{\sigma _{ij} |
| 76 |
> |
}}{{r_{ij} }}} \right)^6 } \right], \\ |
| 77 |
> |
V_{dp} (r_{ij} ,\Omega _i ,\Omega _j ) &= & |
| 78 |
> |
\frac{|\mu_i||\mu_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[ |
| 79 |
> |
\hat{u}_{i} \cdot \hat{u}_{j} - 3(\hat{u}_i \cdot \hat{\mathbf{r}}_{ij}) % |
| 80 |
> |
(\hat{u}_j \cdot \hat{\mathbf{r}}_{ij}) \biggr],\\ |
| 81 |
> |
V_{sticky} (r_{ij} ,\Omega _i ,\Omega _j ) &=& v_0 [s(r_{ij} |
| 82 |
> |
)w(r_{ij} ,\Omega _i ,\Omega _j ) + s'(r_{ij} )w'(r_{ij} ,\Omega _i |
| 83 |
> |
,\Omega _j )].\\ |
| 84 |
> |
\end{eqnarray*} |
| 85 |
|
where $v_0$ is a strength parameter, $s$ and $s'$ are cubic |
| 86 |
< |
switching functions, while $w$ and $w'$ are responsible for the |
| 86 |
> |
switching functions, while $w$ and $w'$ are responsible for the |
| 87 |
|
tetrahedral potential and the short-range correction to the dipolar |
| 88 |
|
interaction respectively. |
| 89 |
|
\[ |
| 92 |
|
w'(r_{ij} ,\Omega _i ,\Omega _j ) = (\cos \theta _{ij} - 0.6)^2 (\cos \theta _{ij} + 0.8)^2 - w_0 \\ |
| 93 |
|
\end{array} |
| 94 |
|
\] |
| 95 |
< |
Although dipole-dipole and sticky interactions are more |
| 95 |
> |
Although the dipole-dipole and sticky interactions are more |
| 96 |
|
mathematically complicated than Coulomb interactions, the number of |
| 97 |
|
pair interactions is reduced dramatically both because the model |
| 98 |
|
only contains a single-point as well as "short range" nature of the |
| 99 |
< |
higher order interaction. |
| 99 |
> |
more expensive interaction. |
| 100 |
|
|
| 101 |
|
\subsection{\label{lipidSection:coarseGrained}The Coarse-Grained Phospholipid Model} |
| 102 |
|
|
| 114 |
|
\caption[A representation of coarse-grained phospholipid model]{A |
| 115 |
|
representation of coarse-grained phospholipid model. The lipid |
| 116 |
|
headgroup consists of $\text{{\sc PO}}_4$ group (yellow), |
| 117 |
< |
$\text{{\sc NC}}_4$ group (blue) and a united $\text{{\sc C}} atom |
| 118 |
< |
(gray) $ with a dipole, while the glycerol backbone includes dipolar |
| 119 |
< |
$\text{{\sc CE}}$ (read) and $\text{{\sc CK}}$ (pink) groups. Alkyl |
| 120 |
< |
groups in hydrocarbon chains are simply represented by gray united |
| 121 |
< |
atoms.} \label{lipidFigure:coarseGrained} |
| 117 |
> |
$\text{{\sc NC}}_4$ group (blue) and a united C atom (gray) with a |
| 118 |
> |
dipole, while the glycerol backbone includes dipolar $\text{{\sc |
| 119 |
> |
CE}}$ (red) and $\text{{\sc CK}}$ (pink) groups. Alkyl groups in |
| 120 |
> |
hydrocarbon chains are simply represented by gray united atoms.} |
| 121 |
> |
\label{lipidFigure:coarseGrained} |
| 122 |
|
\end{figure} |
| 123 |
|
|
| 124 |
|
Accurate and efficient computation of electrostatics is one of the |
| 132 |
|
\] |
| 133 |
|
where $N_A$ and $N_B$ are the number of point charges in the two |
| 134 |
|
molecular species. Originally developed to study ionic crystals, the |
| 135 |
< |
Ewald summation method mathematically transforms this |
| 136 |
< |
straightforward but conditionally convergent summation into two more |
| 137 |
< |
complicated but rapidly convergent sums. One summation is carried |
| 138 |
< |
out in reciprocal space while the other is carried out in real |
| 139 |
< |
space. An alternative approach is a multipole expansion, which is |
| 140 |
< |
based on electrostatic moments, such as charge (monopole), dipole, |
| 144 |
< |
quadruple etc. |
| 135 |
> |
Ewald sum method mathematically transforms this straightforward but |
| 136 |
> |
conditionally convergent summation into two more complicated but |
| 137 |
> |
rapidly convergent sums. One summation is carried out in reciprocal |
| 138 |
> |
space while the other is carried out in real space. An alternative |
| 139 |
> |
approach is the multipole expansion, which is based on electrostatic |
| 140 |
> |
moments, such as charge (monopole), dipole, quadrupole etc. |
| 141 |
|
|
| 142 |
|
Here we consider a linear molecule which consists of two point |
| 143 |
|
charges $\pm q$ located at $ ( \pm \frac{d}{2},0)$. The |
| 163 |
|
$\frac{{d^2 }}{4}$ inside the square root of the denominator is |
| 164 |
|
neglected. This is a reasonable approximation in most cases. |
| 165 |
|
Unfortunately, in our headgroup model, the distance of charge |
| 166 |
< |
separation $d$ (4.63 $\AA$ in PC headgroup) may be comparable to |
| 167 |
< |
$r$. Nevertheless, we could still assume $ \cos \theta \approx 0$ |
| 168 |
< |
in the central region of the headgroup. Using Taylor expansion and |
| 166 |
> |
separation $d$ (4.63 \AA in PC headgroup) may be comparable to $r$. |
| 167 |
> |
Nevertheless, we could still assume $ \cos \theta \approx 0$ in |
| 168 |
> |
the central region of the headgroup. Using Taylor expansion and |
| 169 |
|
associating appropriate terms with electric moments will result in a |
| 170 |
|
"split-dipole" approximation: |
| 171 |
|
\[ |
| 201 |
|
methods. |
| 202 |
|
\begin{figure} |
| 203 |
|
\centering |
| 204 |
< |
\includegraphics[width=3in]{split_dipole.eps} |
| 204 |
> |
\includegraphics[width=\linewidth]{split_dipole.eps} |
| 205 |
|
\caption[Comparison between electrostatic |
| 206 |
|
approximation]{Electrostatic potential map for two pairs of charges |
| 207 |
|
with different alignments: (a) illustration of different alignments; |
| 215 |
|
|
| 216 |
|
\subsection{One Lipid in Sea of Water Molecules} |
| 217 |
|
|
| 218 |
< |
To exclude the inter-headgroup interaction, atomistic models of one |
| 219 |
< |
lipid (DMPC or DLPE) in sea of water molecules (TIP3P) were built |
| 220 |
< |
and studied using atomistic molecular dynamics. The simulation was |
| 221 |
< |
analyzed using a set of radial distribution functions, which give |
| 222 |
< |
the probability of finding a pair of molecular species a distance |
| 223 |
< |
apart, relative to the probability expected for a completely random |
| 224 |
< |
distribution function at the same density. |
| 218 |
> |
To tune our parameters without the inter-headgroup interactions, |
| 219 |
> |
atomistic models of one lipid (DMPC or DLPE) in sea of water |
| 220 |
> |
molecules (TIP3P) were built and studied using atomistic molecular |
| 221 |
> |
dynamics. The simulation was analyzed using a set of radial |
| 222 |
> |
distribution functions, which give the probability of finding a pair |
| 223 |
> |
of molecular species a distance apart, relative to the probability |
| 224 |
> |
expected for a completely random distribution function at the same |
| 225 |
> |
density. |
| 226 |
|
|
| 227 |
|
\begin{equation} |
| 228 |
|
g_{AB} (r) = \frac{1}{{\rho _B }}\frac{1}{{N_A }} < \sum\limits_{i |
| 235 |
|
\end{equation} |
| 236 |
|
|
| 237 |
|
From Fig.~\ref{lipidFigure:PCFAtom}, we can identify the first |
| 238 |
< |
solvation shell (3.5 $\AA$) and the second solvation shell (5.0 |
| 239 |
< |
$\AA$) from both plots. However, the corresponding orientations are |
| 238 |
> |
solvation shell (3.5 \AA) and the second solvation shell (5.0 \AA) |
| 239 |
> |
from both plots. However, the corresponding orientations are |
| 240 |
|
different. In DLPE, water molecules prefer to sit around $\text{{\sc |
| 241 |
|
NH}}_3$ group due to the hydrogen bonding. In contrast, because of |
| 242 |
|
the hydrophobic effect of the $ {\rm{N(CH}}_{\rm{3}} |
| 244 |
|
in DMPC is around the $\text{{\sc PO}}_4$ Group. When the water |
| 245 |
|
molecules are far from the headgroup, the distribution of the two |
| 246 |
|
angles should be uniform. The correlation close to center of the |
| 247 |
< |
headgroup dipole in both plots tell us that in the closely-bound |
| 247 |
> |
headgroup dipole in both plots tells us that in the closely-bound |
| 248 |
|
region, the dipoles of the water molecules are preferentially |
| 249 |
|
anti-aligned with the dipole of headgroup. When the water molecules |
| 250 |
|
are far from the headgroup, the distribution of the two angles |
| 323 |
|
period of 2 ns. During the equilibration period, the system was |
| 324 |
|
initially simulated at constant volume for 1 ns. Once the system was |
| 325 |
|
equilibrated at constant volume, the cell dimensions of the system |
| 326 |
< |
was relaxed by performing under NPT conditions using Nos¨¦-Hoover |
| 326 |
> |
was relaxed by performing under NPT conditions using Nos\'{e}-Hoover |
| 327 |
|
extended system isothermal-isobaric dynamics. After equilibration, |
| 328 |
|
different properties were evaluated over a production run of 5 ns. |
| 329 |
|
|
| 357 |
|
electron density is in the hydrocarbon region. As a good |
| 358 |
|
approximation to the thickness of the bilayer, the headgroup spacing |
| 359 |
|
, is defined as the distance between two peaks in the electron |
| 360 |
< |
density profile, calculated from our simulations to be 34.1 $\AA$. |
| 360 |
> |
density profile, calculated from our simulations to be 34.1 \AA. |
| 361 |
|
This value is close to the x-ray diffraction experimental value 34.4 |
| 362 |
< |
$\AA$\cite{Petrache1998}. |
| 362 |
> |
\AA\cite{Petrache1998}. |
| 363 |
|
|
| 364 |
|
\begin{figure} |
| 365 |
|
\centering |
| 421 |
|
|
| 422 |
|
\section{\label{lipidSection:Conclusion}Conclusion} |
| 423 |
|
|
| 424 |
< |
Atomistic simulations are used in this study to determine the |
| 424 |
> |
Atomistic simulations have been used in this study to determine the |
| 425 |
|
preferred orientation and location of water molecules relative to |
| 426 |
|
the location and orientation of the PC and PE lipid headgroups. |
| 427 |
< |
Based on the result from all-atom simulations, we developed a simple |
| 428 |
< |
coarse-grained model capturing essential features of the headgroup |
| 429 |
< |
solvation which is crucial to transport process in membrane system. |
| 430 |
< |
In addition, the model has been explored in a bilayer system which |
| 431 |
< |
is shown to have reasonable electron density profile, |
| 432 |
< |
$\text{S}_{\text{{\sc cd}}}$ order parameter and other structural |
| 433 |
< |
properties. The accuracy of this model is achieved by matching |
| 434 |
< |
atomistic result. It is also easy to represent different |
| 435 |
< |
phosphorlipids by changing a few parameters of the model. Another |
| 427 |
> |
Based on the results from our all-atom simulations, we developed a |
| 428 |
> |
simple coarse-grained model which captures the essential features of |
| 429 |
> |
the headgroup solvation which is crucial to transport process in |
| 430 |
> |
membrane system. In addition, the model has been explored in a |
| 431 |
> |
bilayer system was shown to have reasonable electron density |
| 432 |
> |
profiles, $\text{S}_{\text{{\sc cd}}}$ order parameter and other |
| 433 |
> |
structural properties. The accuracy of this model is achieved by |
| 434 |
> |
matching atomistic result. It is also easy to represent different |
| 435 |
> |
phospholipids by changing a few parameters of the model. Another |
| 436 |
|
important characteristic of this model distinguishing itself from |
| 437 |
|
other models\cite{Goetz1998,Marrink2004} is the computational speed |
| 438 |
< |
gaining by introducing short range electrostatic approximation. |
| 438 |
> |
gained by introducing a short range electrostatic approximation. |
| 439 |
|
Further studies of this system using z-constraint method could |
| 440 |
|
extend the time length of the simulations to study transport |
| 441 |
|
phenomena in large-scale membrane systems. |
