This report addresses the following problems associated
with the generation of computer models of phospholipid bilayer
membranes using molecular dynamics simulations: arbitrary initial
structures and short equilibration periods, an Ewald-induced strong
coupling of phospholipids, uncertainty regarding which value should be used for surface tension to alleviate the problem of the small size of
the membrane, and simultaneous realization of both order parameters and
the surface area. We generated a computer model of the
liquid-crystalline L-
-dimyristoylphosphatidylcholine
(DMPC) bilayer, starting from a configuration based on a crystal
structure (rather than from an arbitrary structure). To break the
crystalline structure, a 20-ps high-temperature pulse of 510 K (but not
450 or 480 K) was effective. The system finally obtained is an all-atom model, with Ewald summation to evaluate Coulombic interactions and a
constant surface tension of 35 dynes/cm/water-membrane interface, equilibrated for 12 ns (over 50 ns total calculation time), which reproduces all of the experimentally observed parameters examined in
this work. Furthermore, this model shows the presence of significant orientational correlations between neighboring alkyl chains and between
shoulder vectors (which show the orientations of the lipids about their
long axes) of neighboring DMPCs.
 |
INTRODUCTION |
One of the most important structural
characteristics of a biological membrane is that it is like a
two-dimensional liquid; the constituent molecules can undergo thermal
movement within the membrane without altering its overall morphology.
This fluid-like nature of the membrane is largely due to the
conformational freedom of phospholipid molecules and particularly
gauche-trans isomerization of the alkyl chains. There is a
growing body of experimental data on alkyl chain conformational
dynamics in membrane phospholipids.
However, experimental approaches have been insufficient for
understanding many important basic issues, such as the mechanisms by
which small molecules or drugs permeate across the membrane (Subczynski
et al., 1989
, 1990, 1991, 1994; Mason et al., 1991; Wimley and White,
1993; Ashikawa et al., 1994; Peck et al., 1995; Xiang and Anderson,
1995; Paula et al., 1996, 1998; Huster et al., 1997; Frézard and
Garnier-Suillerot, 1998; Mouritsen and Jorgensen, 1998) and the
relationship between the chemical structures of constituent lipid
molecules and the properties of the membrane. For such an
understanding, clearer and more concrete images of the conformations of
lipid molecules in the membrane, i.e., conformations at an atomic
resolution with a time resolution better than the rate of
conformational changes, are required. In addition, with such spatial
and temporal resolutions, new insights may be gained into longstanding
questions in membrane biology, such as the elementary steps of the
translational diffusion of lipid in membranes and the motional
correlation between neighboring lipids and alkyl chains.
One of the most powerful methods for obtaining such information is
molecular dynamics (MD) simulation. Various dynamic processes and
interactions in the membrane have been investigated using MD
simulations, including gauche-trans isomerization; diffusion of small molecules, ions, and water; interaction of lipids with water,
drugs, and peptides/proteins; and the effect of the presence of
cholesterol and unsaturation on phospholipid conformations (Venable et
al., 1993; Pastor, 1994; Merz and Roux, 1996; Jakobsson, 1997; Merz,
1997; Tieleman et al., 1997; Pasenkiewicz-Gierula et al., 2000, and
references therein). However, despite impressive success in these
studies, many problems still need to be addressed regarding the
generation of proper computer models of biological membranes using MD
simulations. These include the use of an arbitrary initial structure
and a short equilibration period, the use of united atoms rather than
all-atom models, a cutoff for Coulombic interactions, and uncertainty
regarding which value should be used for surface tension for a very
small computer model of the membrane. The major objectives in the
present research were to address these problems.
Among these problems, one of the major issues associated with almost
all previous simulations is that the initial structure of the membrane
is produced rather arbitrarily, and an "equilibration" calculation
is performed for only a very short time, in most cases for less than a
few nanoseconds. This length of time is so short that neither the alkyl
chain conformation nor the orientation of a lipid about its long axis
has time to reach an equilibrium value (as will be shown in this
report, the relaxation times for these processes are substantially
greater than 2 ns). Therefore, the structure of the membrane model even
at the time of the production run is, in many cases, dominated by the
arbitrary initial structure. For example, in most simulations, the
orientation of the phospholipid about its long axis is randomized in
the initial structure, but, as shown in this paper, the orientation is
actually not totally random; neighboring phospholipids tend to show
correlated orientations due to the presence of two alkyl chains in a
phospholipid molecule, which is difficult to reproduce in a simulation
without prolonged (at least several nanoseconds of) equilibration periods.
In the present study, we started MD calculations with an initial
structure based on a crystal structure and followed the development of
the system toward equilibrium for 13 ns. The total calculation time
spent for various trials of simulation conditions exceeded 50 ns. We
often found that several nanoseconds were required for the system to
reach equilibrium after the simulation conditions, such as temperature
and surface tension, were changed.
Initially, a united-atom approach and an 8-Å cutoff for both van der
Waals and Coulomb interactions were used (system 1). After the system
was judged to have reached equilibrium under these conditions (5.2 ns
after beginning the calculation), a production run was carried out for
800 ps. At 6.0 ns, the cutoff for Coulombic interactions was abolished
and an Ewald method was introduced for a more precise evaluation of
Coulombic interactions. At 8.2 ns, surface tension was included, and
the united-atom model was changed to an all-atom model (system 2). The
final system was composed of 112 dimyristoylphosphatidylcholine (DMPC)
molecules and 3016 water molecules (total of 22,264 atoms) and was
found to be stable under constant temperature, pressure, and surface tension.
Previously, Tu et al. (1995) approached similar problems and carried
out a 1550-ps constant pressure/temperature MD simulation of a bilayer
of dipalmitoylphosphatidylcholine with an all-atom model, starting from
the crystal structure of phosphatidylcholine and using the all-atom
potentials they developed. They showed that the bilayer was stable
throughout the 1550-ps simulation and that it well reproduced many
important experimental observations. In the present study, we used an
AMBER all-atom force field (Cornell et al., 1995), with an extended
ensemble, including surface tension, which was applied to circumvent
the problems that might be caused by the small size of the simulation
system. We carried out an extensive examination of various simulation
conditions. In the present report, we concentrate on describing the
development of computer models of the DMPC membrane. A detailed
analysis using the obtained membrane will be presented elsewhere.
| |
METHODS |
Two fully hydrated DMPC membranes (system 1 and system 2) were
generated. We initially generated system 1, a united-atom model, from a
crystal structure, using a cutoff distance of 8 Å (for both Coulomb
and van der Waals interactions), but not surface tension. System 2 was
constructed from the fully equilibrated system 1. It used an all-atom
model, Ewald summation to evaluate Coulombic interactions, a 15-Å
cutoff distance for van der Waals interactions, and a 70-dyn/cm surface
tension for the bilayer. These simulation conditions are
summarized as a function of time in Fig.
1.
View larger version (19K):
[in this window]
[in a new window]
|
FIGURE 1
The time schedule of simulation conditions and
temperature variations used in this simulation. During the first 6.0 ns, nonbonded interactions were evaluated using a cutoff at 8 Å. After
6.0 ns, Coulomb interactions were evaluated correctly using an Ewald
summation, and van der Waals interactions were cut off at 15 Å.
Surface tension (Zhang et al., 1995) was applied to system 1 at 9.5 ns.
Meanwhile, system 2, an all-atom model, was constructed at 8.2 ns and
included Ewald evaluation of the Coulombic interactions and the surface
tension. Eight test runs were carried out for each system: after 9.5 ns
for system 1 and after 8.2 ns for system 2, totaling over 35 ns. Time
for system 2 was counted separately, starting at 8.2 ns (which is set
at 0 ns).
|
|
Initial structure of system 1
L--Dimyristoylphosphatidylcholine (DMPC; Fig.
2) was constructed based on the atomic
coordinates of 1,2-dilauroylphosphatidylethanolamine (DLPE) (Hitchcock
et al., 1974). First, to generate myristoyl chains, two
-CH2- groups in a trans conformation
were inserted before the terminal methyl groups in the lauroyl chains.
Second, to produce the choline group, the hydrogen atoms of the amine group of DLPE were replaced with the methyl groups. The structure was
then optimized using AMBER 4.0 software (Pearlman et al., 1991). To
reduce the computation time, a united-atom approximation, in which -CH,
-CH2, and -CH3 groups were
treated as single united atoms, was used for system 1. Therefore, each
DMPC molecule consisted of 46 particles.
View larger version (20K):
[in this window]
[in a new window]
|
FIGURE 2
Numbering of atoms and torsion angles for the DMPC
molecule according to Sundaralingam (1972). Atoms other than P, O, and
N are carbon atoms with 1, 2, or 3 hydrogen atoms attached.
|
|
The DMPC bilayer was constructed by arranging 56 DMPC molecules
according to the symmetry of the crystal structure of DMPC (Pearson and
Pascher, 1979; no atomic coordinates were given in this report). In
each leaflet, 28 DMPC molecules were placed in a two-dimensional 4×7
array (see Fig. 4, initial). The bilayer was immersed in a
64-Å × 48-Å × 38.5-Å box filled with TIP3P water (Jorgensen et
al., 1983). Water molecules whose oxygen and hydrogen atoms were closer
than 3.0 Å and 2.15 Å, respectively, to any DMPC atom were removed.
The total number of water molecules in the system was 1197, and the
total number of atoms was 6167. Approximately 10-Å-thick water layers
were formed on both sides of the lipid bilayer, with 80 water molecules
left in the hydrophobic region of the membrane. The location of water
molecules was then optimized using AMBER (Pearlman et al., 1991) with
the position of DMPC molecules constrained. The final dimension of the
system was 57.5 Å × 45.7 Å × 37.5 Å.
Initial structure of system 2
The initial structure of system 2 was constructed based on the
structure of fully equilibrated system 1 after 8.2 ns of MD simulation
(Fig. 1). The united-atom approximation, used for system 1, was changed
to an all-atom model. To circumvent problems that might be caused by
the small size of the system (56 DMPCs), the number of DMPC molecules
in system 2 was doubled. This was done by placing two sets of the
equilibrated system 1 (47.97 Å × 32.24 Å, 61.19 Å thick) next to
each other. The final dimensions were 47.97 Å × 64.48 Å. In
addition, a 3-Å-thick water layer was added to each side of the
membrane to further ensure full hydration. In system 2, the number of
DMPC molecules was 112, and the number of water molecules was 3016 (1197 × 2 + 622). The total number of atoms was 22,264.
Atomic parameters
The charge distribution on the DMPC molecule was calculated
using the GAUSSIAN-92 program (Frisch et al., 1992). For system 1, after ab initio calculation with the STO-3G basis set, atomic charges were determined by fitting to the electrostatic potential calculated at the points selected according to the Merz-Singh-Kollman scheme (Singh and Kollman, 1984; Besler et al., 1990) (Table
1). For
system 2, the 6-31G* basis set and restrained electrostatic potential
(RESP) fitting (Bayly et al., 1993) were used. To calculate atomic
charges for system 2, four representative conformations of lipid
molecules in system 1 were selected, based on the distances between the
nitrogen atom in the headgroup and two glycerol oxygens (O22 and O32,
see Fig. 2, Stouch and Williams, 1992; Pasenkiewicz-Gierula et al.,
1997). Nitrogen atoms are widely spread along the axis normal to the
bilayer and since the distributions of the distances between the
headgroup nitrogen and two glycerol oxygens showed four peaks (Stouch
and Williams, 1992; Pasenkiewicz-Gierula et al., 1997), a
representative conformation of DMPC was selected from each peak. We
carried out ab initio calculations for the four different conformations
of DMPC and then calculated the atomic charge set to reproduce four
sets of electrostatic potentials (Stouch and Williams, 1992;
Pasenkiewicz-Gierula et al., 1997).
To calculate other interactions (bond stretching, bending, torsion, and
van der Waals), the OPLS parameter set (Jorgensen and Tirado-Rives,
1988) was used for system 1 and the AMBER parameter set (Cornell et
al., 1995) was used for system 2.
Simulation conditions
The MD simulation was carried out using AMBER 4.0 (Pearlman et
al., 1991) and AMBER 4.1 (Pearlman et al., 1995) for systems 1 and 2, respectively. In-house subroutines for Ewald summation and the
Zhang-Nosé algorithm (see below) were added to the AMBER 4.1 package. The SHAKE algorithm (Ryckaert et al., 1977) was used. The time
step was set at 2 fs. Under a periodic boundary condition, the
temperature and pressure were kept constant during the simulation of
system 1 (Berendsen et al., 1984; Nosé, 1984; Zhang et al., 1995).
At the initial stages of the simulation, a residue-based cutoff was
used for nonbonded interactions with a cutoff distance of 8 Å. DMPC
molecules were considered to have three residues, i.e., a
phosphorylcholine plus glycerol group and two alkyl chains. After the
DMPC conformation had reached an equilibrium state (at 6.0 ns, see Fig.
1 and Results and Discussion for details), the cutoff algorithm for
Coulombic interactions was changed to an Ewald method (Ewald, 1921).
The convergence parameter was set at 0.150 Å
1. Integer vector K
values in the reciprocal space were taken to K 2 = 314 (the average value during the period
when the trajectory of the equilibrated system 2 was obtained (1.0-1.5
ns after the system 2 simulation started)), which varied according to
changes in the cell size. Van der Waals interactions were truncated at 15 Å.
The temperature of the system was elevated first to break the
crystalline structure of the bilayer and then to accelerate the
approach to the equilibrium state. The scheme of the temperature variations used in this study is also shown in Fig. 1. Finally, the
temperature was set at 310 K and the system was equilibrated.
Application of the Ewald method to system 1 (6-13 ns) caused a
decrease in the surface area per lipid and the lateral movement of
lipid molecules (by 9.5 ns). To counteract these changes, surface tension at the membrane-water interface (Zhang et al., 1995) was instituted at 9.5 ns. The Zhang method (1995) was combined with the
Nosé method (1984) to control the temperature, which enabled us
to carry out MD simulation in an NP
T extended ensemble, where denotes the surface tension. In this scheme, the Hamiltonian of the
system is given as
|
(1)
|
where pi and
ps are the momentum of particle
i and the additional degree of freedom s,
respectively; mi is the mass of
particle i; Q and W are the mass for
the motion of s and volume motion, respectively;
g is the number of degrees of freedom; h is
the cell size; V is the volume; and
T0,
P0, and 0 are
the reference temperature, pressure, and surface tension, respectively.
Various values for surface tension (0, 3, 5, 10, 15, 20, 25, 30, 35, and 50 dyn/cm per membrane-water interface) were tested to find one for
which both surface area per lipid and the order parameter profile
agreed with experimental values. For many values of the surface
tension, the approach to equilibrium after a change in the surface
tension value required a few nanoseconds. However, no proper value for
the surface tension could be found that could satisfy both the surface
area/lipid and the order parameter profiles. We assumed that this could
be caused by the united-atom approximation and/or the size of the
system. Therefore, we constructed a new DMPC bilayer model system 2 containing twice as many DMPC molecules with the all-atom model rather
than the united-atom model, still using surface tension and Ewald
summation for Coulombic interactions and a 15-Å cutoff for van der
Waals interactions. After examining various values for surface tension
in system 2, a value of 35 dyn/cm/interface (70 dyn/cm/bilayer) was
found to generate a stable membrane and satisfied our criterion of
satisfying both the surface area/lipid and the order parameter profile
across the bilayer.
MD simulation of system 1 was carried out on a SPARCstation 10 with a
hardware accelerator for molecular dynamics simulations (MD Engine;
Amisaki et al., 1995; Kitamura et al., 1995; Toyoda et al., 1995,
1999). For calculations of system 2, Sparc CPU-5CE Force
computers were used.
| |
RESULTS AND DISCUSSION |
This section consists of two subsections. We first describe the
process of system development from the crystalline structure to
equilibrium (section 1). Initially, we used a united-atom model and an
8-Å cutoff for van der Waals and Coulomb interactions (system 1).
After equilibration, based on the equilibrated structure of system 1, we constructed the next membrane system (system 2), in which the size
of the membrane was doubled, and an all-atom model, an Ewald approach
to evaluate Coulomb interactions, a constant surface tension, and a
15-Å cutoff for van der Waals interactions were used. The reasons we
wished to develop system 2 and the details of the development and
approach toward equilibrium of system 2 are described in the latter
part of section 1. In section 2, we characterize the fully equilibrated
system 2, in comparison with the results of experiments and other simulations.
Equilibration of system 1
Definition of three parameters that were used to characterize the
membrane
To examine the approach of system 1 to the thermally equilibrated
structure, we mainly monitored three parameters which represent (1) the
order of the alkyl chains, (2) reorientation about the long axis of the
molecule, and (3) the relative positions of the lipid molecules in the
membrane plane. In the following discussion, 
ensemble denotes an ensemble average. Ensemble
averaging at a given time point involves averaging over all lipid
molecules and over a 50-ps time interval (50 consecutive data points 1 ps apart, from (t 
49) ps until t). This
time averaging over 50 ps is useful to suppress fast noise but has
nothing to do with ensemble averaging over the phase space.
The order parameter for a given segment of an alkyl chain was estimated
as follows. An instantaneous orientation of the alkyl chain at the
ith carbon is defined by a segment vector
Ci1 
Ci+1, which is
designated as ci. An instantaneous angle
i is defined as the angle between
ci and
ciensemble. An order parameter
Smol(i) can be expressed as
(Hubbell and McConnell, 1971)
|
(2)
|
As the average tilt angle of the long axes of alkyl chains is
nearly 0° in the liquid-crystalline membrane (Meier et al., 1982;
Moser et al., 1989; reconfirmed in this work as described later),
ciensemble is almost parallel
to the membrane normal after equilibration. Hence, for determination of
Smol in the liquid crystalline
membrane after equilibration, we employed as i
the angle between ci and
ciensemble the membrane
normal, which agrees with the conventional definition of
Smol in the NMR literature (Seelig and
Niederberger, 1974). During the equilibration period when the average
tilt angle of the alkyl chains are large, the definition of
i as the angle between
ci and was employed, and the order
parameter using this definition is expressed as 
.
This definition allows independent evaluation of the cooperative tilt
of alkyl chains and the instantaneous deviation of each segment's
orientation from the tilt axis, both of which are important parameters
to assess the system's approach toward equilibrium. To describe the
approach to equilibrium, we also use an averaged order parameter of a
chain, which is defined as the mean value of the order parameter values
() over 10 segments (from carbon numbers 4 to 13).
Alkyl chain orientation t is defined as the unit vector
pointing from the midpoint of C1 and C2 toward that of C13 and C14, and
alkyl chain tilt is defined as the angle between t and the
membrane normal.
According to Cornell et al. (1995), the alkyl chain torsion angle
potential has three minimum values at 68°, 180°, and 68° (292°). These are generally noted as gauche plus
(g+), trans (t), and
gauche minus (g),
respectively. Actual angular distributions of several torsion angles
observed in system 2 are shown in Fig. 3,
which shows that the distribution for both trans and
gauche conformations are broad. To cover the entire range of
the distribution, in the present report, the torsion angles are
classified as g+, t, and
g when they are between 0° and
120°, 120° and 240°, and 240° and 360°, respectively. Kink refers to a conformation in which the three adjacent torsion angles are
g+tg or
gtg+.
View larger version (25K):
[in this window]
[in a new window]
|
FIGURE 3
Angular distributions of torsion angles observed in
system 2. Dotted line,  4; thin solid line, 9; thick solid line,
14. See Fig. 2 for the definitions of the numbers of torsion
angles.
|
|
The second parameter reflects the distribution of the reorientation
angle about the molecular long axis. A vector connecting O21 and C3 of
each DMPC (see Fig. 2, shoulder vector) is projected on the membrane
surface, and the angle made by this vector at time t and
time 0 (initial arrangement) is defined as the reorientation angle
(t). For each time point, (t) was averaged
over the last 50-ps interval.
The third parameter monitors the relative positions of DMPC molecules
in the membrane plane, using a two-dimensional Delaunay tessellation
analysis (Takaoka and Kitamura, 1996; Shinoda and Okazaki, 1998). In
this analysis, closest DMPC neighbors using the mass centers of lipid
molecule projected on the membrane plane are connected by straight
lines, which form Delaunay triangles. These triangles cover the entire
membrane plane, and their arrangement has been shown to be unique
(Tanemura et al., 1983). This approach has been applied in the analysis
of amorphous systems and protein conformations (Yamato et al., 1994;
Gerstein et al., 1995), and is suitable for describing changes in the
local arrangement and topography of particles.
In system 1, 56 triangles were formed for each monolayer consisting of
28 DMPC molecules. As DMPC molecules diffuse in the membrane plane,
their relative positions and nearest neighbors change, and so does the
local arrangement of the triangles. As each triangle is identified by
the three molecules at the vertices, the extent of molecular
rearrangement can be traced by counting the number of triangles that
represent the original set of three lipids.
We first constructed this triangle set at 0.85 ps and then
reconstructed it every 100 ps. By counting the number of triangles that
remained with unchanged lipids, we monitored the degree of DMPC
rearrangement. When an original (initial) triangle was reformed after
it had disappeared, it was counted as one of the remaining triangles
even if it had been broken for a long period. The major cause for the
formation of new triangles is not the exchange of the lipids'
locations (i.e., translational diffusion) but small readjustment of
their locations, i.e., local sliding around of the lipids, which
reflects local expansion and compression of the membrane.
Equilibration stage 1 (0-650 ps): water exclusion from the
membrane
Immediately after initiation of this simulation, large
translational motion of DMPC molecules in the membrane plane took
place, and the chain conformation started to lose its initial
structure. The bilayer structure remained stable. After 40 ps of
equilibration, the temperature was set at 310 K. The scheme for the
temperature setting over time used in this simulation is shown in Fig.
1.
The water molecules that had been present in the hydrophobic region of
the membrane in the initial structure showed remarkable movements
during this period (Fig. 4). First, they
assembled in the central part of the bilayer and formed two clusters
during the first 20 ps, and then they formed linear chains parallel
to the alkyl chains, in which water molecules were bonded via hydrogen bonds (Fig. 4, 20-50 ps). These water molecules left the hydrocarbon region of the bilayer one by one in a line (Fig. 4, 100-50 ps). Thereafter, none of the water molecules reached deep inside the membrane during this simulation (total time over 50 ns). The
experimentally observed permeability of water across the membrane is
large considering the hydrophilic nature of water, and its mechanism is
not yet understood (Subczynski et al., 1994). The result obtained here suggests that water molecules move across the membrane as clusters or
lines rather than as single molecules, which would dramatically increase the water permeation rate. Such passage may be enhanced in the
presence of large cooperative movement of lipids in the membrane.
View larger version (88K):
[in this window]
[in a new window]
|
FIGURE 4
Snapshots of the system during its early equilibration
stages, showing that water molecules that had been inside the membrane
in the initial structure moved out of the hydrophobic region of the
membrane. DMPC molecules are drawn in thin lines, while water is shown
using a space-filling model. Red denotes oxygen atom, blue nitrogen,
yellow phosphate, white hydrogen, and green united carbon atoms. Water
molecules first assembled in the central part of the bilayer to form
two clusters by 20 ps and then formed lines parallel to the alkyl
chains, in which water molecules were bonded via hydrogen bonding
(20-50 ps). They went out of the bilayer in a line (100-150
ps).
|
|
Stage 2 (650 ps to 1.2 ns): high temperature pulses abolish alkyl
chain tilt
In the initial structure, alkyl chains were tilted ~30°C with
respect to the bilayer normal (Fig. 4). This tilt persisted even after
650 ps. Experimental results have indicated that the average tilt angle
of alkyl chains in DMPC bilayers is nearly 0°C in the
liquid-crystalline phase (Meier et al., 1982; Moser et al., 1989). It
was possible that our system could not reach a state of thermal
equilibrium at 310 K during the available calculation time due to the
high activation energy required to change the collective tilt angle in
the crystal structure. This difficulty could be enhanced by the use of
a periodic boundary condition (too much coupling between DMPC molecules).
Therefore, the system temperature was raised to 450 K for 20 ps (660 ps ~ 680 ps), which was considered to be sufficient to break the
crystalline structure of the membrane, and then slowly lowered to 310 K. However, this maneuver did not induce large changes in the system.
As shown in Table 2, both the tilt angle and the average order parameter of the alkyl chains were still much
greater than the experimental values.
View this table:
[in this window]
[in a new window]
|
TABLE 2
Time development of the mean tilt angle of alkyl chains
with respect to the membrane normal, average order parameter, and the
number of gauche and kink conformations
per myristoyl chain
|
|
Next, the system temperature was raised again to 510 K for 20 ps
(800-820 ps) and then slowly lowered to 310 K. This effectively changed the mean tilt of the hydrocarbon chains to nearly 0°C and
reduced the order parameters close to the experimental values (Table
2). In a separate test simulation, we found that a temperature pulse of
480 K for 20 ps did not promote such changes. It is somewhat difficult
to explain why such changes occurred at 510 K and not at 480 K or 450 K, but the marked contrast in response to high temperature pulses
between 510 K and 480 K (or 450 K) is in agreement with the fact that
the alkyl chain tilt angle in a bilayer made of a single lipid species
is determined cooperatively.
Stage 3 (1.2-2.5 ns): equilibration in terms of the alkyl chain
order
Fig. 5 shows the order parameter
(mol) at the positions of several carbon
atoms in the -chain plotted as a function of time along with the
temperature profile. Because initially all of the alkyl chains were in
the fully extended (all-trans) conformation, the order
parameter was initially large and then gradually decreased. At 310 K,
between 1.2 and 1.6 ns, the order parameters were still higher than the
corresponding experimental values (Fig. 5). To accelerate the
equilibration process and reduce Smol
values, the temperature of the system was raised to 360 K for 300 ps
(1.6-1.9 ns) and then lowered to 310 K. However, even after this
operation, the order parameters were still considerably greater than
the experimental values. Therefore, the temperature was again raised to
360 K (2.1-4.7 ns) to accelerate equilibration. After 0.4 ns at 360 K
(2.5 ns after beginning the simulation), all of the order parameters
(and thus the average order parameter) reached stable values. (When the
temperature was lowered to 310 K at 4.7 ns and the system was further
equilibrated, the observed profile of the order parameter across
the membrane agreed with the experimental profile.)
View larger version (40K):
[in this window]
[in a new window]
|
FIGURE 5
Order parameters (mol) at
selected carbon atoms in the chain plotted as a function of time.
Solid line, C23; dashed line, C26; dotted line, C29; dot-dashed line,
C212. See text for a definition of the segmental order parameter.
|
|
Stage 4 (2.5-4.7 ns): equilibration regarding lipid reorientation
about its long axis and lipid exchanges
Fig. 6 a shows the
projection of the shoulder vector O21 C3 on the membrane plane for
each DMPC molecule at its mass center at different time points in
the membrane. At 200 ps, the initial 4 × 7 arrangement of DMPC
molecules can be seen clearly, and the shoulder vectors show only small
reorientation. Fig. 6 b shows the angular distribution of
the shoulder vectors with respect to their initial orientation. The
distribution initially centered around 0° and gradually spread to
become flat at ~4.5 ns in the MD simulation. This suggests that no
long-range order exists with regard to the orientation of the shoulder
vector.
View larger version (35K):
[in this window]
[in a new window]
|
FIGURE 6
(a) The vectors connecting
O21C3 projected on the membrane surface (shoulder
vector), showing the reorientation of DMPC molecules about their long
axes. Vectors of 28 molecules in one layer are shown.
(b) Distributions of the DMPC reorientation angle about
the molecular long axis at 200 ps (thin solid line), 600 ps (dashed line), 1.5 ns (dot-dashed
line), 3.0 ns (dotted line), and 4.5 ns
(thick solid line). The shoulder vector for each DMPC
molecule is projected on the membrane surface, and the angle of the
vectors at times t and 0 (initial arrangement) is
defined as the reorientation angle.
|
|
However, whether or not there remain short-range correlations in the
orientation of the shoulder vectors of nearby lipids even after the
system has reached equilibration is a separate issue, which will be
addressed later. The two-dimensional reorientational autocorrelation
function of the shoulder vector (P1(cos)) was adequately fitted to a single exponential function, and the relaxation time of reorientation about the DMPC long axis was estimated to be
~1.5 ns at 360 K. This result suggests that, for full equilibration of the system regarding reorientation about the molecular long axis,
~4.5 ns (three times the relaxation time) will be required, as
confirmed later in the present paper.
Fig. 7 shows the time-dependent decrease
in the number of Delaunay triangles remaining from the initial set (at
0.85 ps). This result shows that ~4.7 ns is required for the local
rearrangement of nearest neighbors. Note that, as pointed out in the
paragraphs describing the Delauney tessellation analysis, this time
does not reflect the frequency of exchange of lipid locations but is rather a yardstick for the time required for relaxation of lipid local
densities. In fact, the trajectories of the mass centers of DMPC
molecules indicate that lipid hop movement is infrequent. The diffusion
constant evaluated from the mean-square displacement of DMPC with
respect to the mass center of the whole system and that with respect to
their nearest neighbors at t = 0 were 4.3 × 107
cm2/s and 3.0 × 107
cm2/s, respectively, which are calculated using a
trajectory of system2 (0.5-5 ns). These values are greater than
experimentally observed diffusion rates by a factor of 3-10 (Vaz et
al., 1985), suggesting that these diffusion constants represent local
jittering motion rather than true translational diffusion, consistent
with the interpretation of Delaunay tessellation analysis.
View larger version (21K):
[in this window]
[in a new window]
|
FIGURE 7
The number of Delaunay triangles remaining at time
t from the number at time 0.85 ps. Delaunay tessellation
(triangle set) is reconstructed every 100 ps.
|
|
Stage 5 (4.7-6.0 ns): equilibrated membrane with a united-atom
approximation and an 8-Å cutoff
The changes in the three observable parameters, the average order
parameter, angular distribution of the shoulder vector, and the number
of remaining Delaunay triangles, over the simulation time are
summarized in Fig. 8. The approximate
times at which these parameters were judged to reach equilibrium values
are indicated. Angular distribution was parameterized by
|
(3)
|
where p() is the distribution of the angle between
the projection of the shoulder vector on the membrane plane and its
initial orientation. R(t) is 1 when
p() becomes constant (homogeneous distribution), and 0 at
time 0.
View larger version (32K):
[in this window]
[in a new window]
|
FIGURE 8
Time course of the changes in the three parameters that
were used to examine the approach of the system to the equilibrated
state (bottom), along with the time schedule of
temperature variations (top). Arrows at the bottom
indicate the times at which these and several other parameters were
judged to reach their respective equilibrium values: solid line,
reorientation angle about the long axis R(t), as defined
in Eq. 2; dashed line, the order parameter ()
averaged over 10 segments (from carbon numbers 4 to 13) in the
-chain; dotted line, the fraction of remaining Delaunay triangles,
representing lipid exchange reactions. The number at time 0.85 ps is
normalized to 1.
|
|
Based on these observations, we concluded that at 360 K the system lost
its memory of the initial arrangement and reached a thermal equilibrium
state at ~4.7 ns of MD simulation (Fig. 8). At this time, the system
temperature was lowered to 310 K and the MD simulation was continued
for an additional 0.8 ns for equilibration at 310 K until 5.5 ns
and then for a production run between 5.5 and 6.0 ns. We used the
last 500-ps trajectory for the characterization of this model (system
1), which will be described later together with system 2. In this
membrane, cooperative tilt of alkyl chains was completely eliminated,
and both the average order parameter and the numbers of
gauche and kink conformations are in the range of
experimental observations (in Table 2, more precise comparison with
experiments is somewhat difficult due to great variations in
experimental results, particularly the data on the number of
kink conformations; see below (Chain Conformation) for details).
Stage 6 (6.0-13.0 ns): elimination of the cutoff for Coulombic
interactions
At 6.0 ns, we stopped using an 8-Å cutoff for nonbonded
interactions and switched to a more accurate calculation, in which Coulombic interactions were evaluated using the Ewald method
(Ewald, 1921, Fig. 1). In addition, the cutoff distance for van der
Waals interactions was increased from 8 Å to 15 Å. We applied the
Ewald method because the use of a cutoff for Coulombic interactions had
been reported to induce various problems in MD simulations (Loncharich
and Brooks, 1989; Kitchen et al., 1990; Schreiber and Steinhauser,
1992; Guenot and Kollman, 1993; Saito, 1994; Oda et al., 1995).
After switching to the Ewald method, we noticed that the lateral
movement of DMPC molecules was reduced and the simulation box began to
shrink in the direction of the membrane plane. Similar phenomena were
reported by Egberts et al. (1994), who noted that their system
approached a gel-like state. This may be due to the application of the
Ewald method to a very small system and can be circumvented by
including the surface tension, which may be useful for simulating a
small membrane system. Therefore, the surface tension was incorporated,
as described in Methods, at 9.5 ns.
Various values for surface tension were examined (0, 3, 5, 10, 15, 20, 25, 30, 35, and 50 dyn/cm per water-lipid interface). An MD simulation
was carried out for 200-4000 ps for each value. We tried to find a
proper value of the surface tension for the present system based on the
criterion that the model membrane should simultaneously reproduce the
D-NMR's Smol profile (Seelig and
Seelig, 1974) and the surface area per DMPC as estimated by NMR and
x-ray diffraction (from 56 to 72 Å2 for DMPC and
DPPC membranes; Büldt et al., 1979; Nagle and Weiner, 1988; Rand
and Parsegian, 1989; Nagle, 1993; Nagle et al., 1996; Koenig et al.,
1997; Petrache et al., 1998; Fig. 9).
However, none of these values for the surface tension satisfied the
above criterion. For example, with 5 dyn/cm/interface, the order
parameter profile was well reproduced, but the surface area/DMPC was
only 57 Å2, which is close to the smallest value
in the distribution of experimental values (56-72
Å2); i.e., this membrane reproduced almost all
of the experimentally obtained parameters, including the order
parameter profile of the myristoyl chain and the relative amounts of
gauche and kink conformations, except for the
surface area/DMPC.
View larger version (24K):
[in this window]
[in a new window]
|
FIGURE 9
Comparison of the surface area per lipid molecule
obtained by several simulations and experiments. Open symbols indicate
the values obtained in this study (m, system 1;
n, system 2). (a) Chiu et al., 1995;
(b) Pasenkiewicz-Gierula et al., 1997;
(c) Berger et al., 1997; (d) Shinoda et
al., 1997; (e) Marrink et al., 1993; (f)
Egberts et al., 1994; (g) Tieleman and Berendsen, 1996;
(h) Marrink et al., 1998; (i) Tu et al.,
1995; (j) Shinoda et al., 1995; (k)
Shinoda and Okazaki, 1998; (l) Smondyrev and Berkowitz,
1999; (o) Büldt et al., 1979; (p)
Rand and Parsegian, 1989; (q) Nagle and Weiner, 1988;
(r) Nagle et al., 1996; (s) Rand and
Parsegian, 1989; (t) Petrache et al., 1998;
(u) De Young and Dill, 1988; (v) Thurmond
et al., 1991; (w) Nagle, 1993.
|
|
With a greater surface tension, the membrane was initially stretched
too much, which first induced alkyl chain tilt and then great decreases
in the order parameters. In some cases, the membrane became unstable.
Without surface tension, the order parameters became too large and the
surface area per DMPC became too small. This sensitivity of the surface
area to surface tension was also recently reported by Feller and Pastor
(1999).
Stage 7 (8.2-12.0 ns or 0-5.4 ns after an all-atom calculation
was initiated): system 2, an all-atom model
The result described above, i.e., realization of the proper order
parameter profile (in fact, proper values for almost all experimental
values), but with less surface area compared with the experimental
value, suggested that this problem may be due to the use of the
united-atom approximation. In this model, the volume and the geometries
of -CH, -CH2, and -CH3
groups cannot be reproduced well, yet these groups represent the
majority in our membrane model. Therefore, we generated an all-atom
model of the DMPC bilayer (system 2) based on the structure of the
fully equilibrated system 1 at 8.2 ns.
After 50 ps (hereafter, time 0 refers to the initial time of system 2 simulation, i.e., 8.2 ns at system 1) of equilibration of system 2 using an NVT ensemble, which is followed by a 200-ps NPT ensemble
simulation, a minus pressure was added parallel to the membrane plane
direction, to stretch the surface area from 58 Å2 to 65 Å2/DMPC, which
is around the middle value of experimental evaluations. An NVT ensemble
simulation was carried out again, until the potential energy and cell
dimension normal to the membrane plane converged to the equilibrium.
At 500 ps, we switched to an extended NPT ensemble to test various
values of surface tension (0, 10, 15, 20, 25, 28, 30, 35, and 40 dyn/cm
per interface) and found that 35 dyn/cm per lipid-water interface
(70 dyn/cm per bilayer) was appropriate for this system. At a surface
tension of 35 dyn/cm, the system appeared stable over 4.9 ns (5.4 ns
for all-atom calculation to 0.5 ns for equilibration calculation before
an extended NPT ensemble was introduced). The profile of the order
parameter across the bilayer was close to that found in deuterium NMR
(see below), and the surface area per DMPC was 61.5 ± 0.6 Å2, which is close to the experimental
observations (Büldt et al., 1979; Nagle and Weiner, 1988; Rand
and Parsegian, 1989; Nagle, 1993; Nagle et al., 1996; Koenig et al.,
1997; Petrache et al., 1998). Thus, we concluded that a surface tension
of 35 dyn/cm/interface (70 dyn/cm/bilayer) satisfies our criteria.
Therefore, a membrane with a surface tension of 70 dyn/cm was used for
further characterization of the membrane.
Feller et al. (1995) calculated the surface tension of a lipid/water
interface based on an MD simulation of a DPPC bilayer membrane at a
constant surface area per lipid (65.1 Å2) and
reported that it increased from 6.3 dyn/cm with a cutoff method to 33.5 dyn/cm with an Ewald summation (Feller and Pastor, 1996). This result
agrees well with the present simulations.
A question then arises as to why, as reported elsewhere (Pastor, 1994;
Merz, 1997; Pasenkiewicz-Gierula et al., 1997, 1999), the united-atom
approach coupled with a cutoff of nonbonded interactions works well in
terms of agreement with experimentally observable parameters and why
inclusion of Ewald summation (and surface tension) in the united-atom
approach does not work well. We propose that the application of Ewald
summation to a very small system under a periodic boundary condition,
such as in the present case, induces very strong coupling between
molecules and that the use of a cutoff helps alleviate this problem. In
addition, it is possible that the OPLS parameters we used in this
simulation may not be suitable for MD simulations with Ewald summation
because they have been optimized for simulations with cutoff of
nonbonded interactions, and they may not be readily useful for MD
simulations with Ewald summation (Andrea et al., 1983; Kitchen et al.,
1990; Smith and Pettitt, 1991; Oda et al., 1995; Hünenberger and
McCammon, 1999).
The total calculation time spent to examine various simulation
conditions (including both the united-atom and all-atom models) exceeded 50 ns. This was necessary for examining simulation parameters and equilibrating the system. During this whole period, several occasions of actual hop exchange of DMPC molecules, which must be an
elementary step of the translational diffusion of lipids, were noted.
However, no flip movement of lipid across the membrane occurred.
Characterization of the membrane (system 2) in the equilibrium
state
Fig. 10 is a snapshot of a fully
equilibrated membrane (system 2) at 3.5 ns. In the present section, we
report analyses carried out using a 1000-ps trajectory of system 2 between 2.9 and 3.9 ns. As atomic coordinates were saved every 1 ps, we
used 1000 time points for analyses. The mean values of the various
parameters reported here represent averages over both time (1000 ps,
1000 points) and all of the molecules (112 DMPCs), except as
specifically stated. For comparison, we present some analyses for
system 1, which were carried out using a 500-ps trajectory between 5.5 and 6.0 ns.
View larger version (135K):
[in this window]
[in a new window]
|
FIGURE 10
Snapshot of the DMPC bilayer-water system at 3.5 ns
(system 2). Red represents oxygen, blue nitrogen, yellow phosphorus,
white hydrogen, and green carbon atoms. Hydrogen atoms attached to DMPC
molecules are omitted for clarity.
|
|
Surface area
The surface area/DMPC of system 2 was 61.5 ± 0.59 Å2 (that of system 1 was 59.3 ± 0.89 Å2), which is in agreement with experimental
results, considering variations of the data (Fig. 9). Feller et al.
(1995) carried out a constant surface area simulation and monitored the
free energy of the system. They concluded that 68.1 Å2 (DPPC) gave the minimum free energy among
four surface areas tested, and the average surface tension of that
system was 74.0 dyn/cm per bilayer (37.0 dyn/cm per interface). Feller
and Pastor (1999) recently applied various values of surface tension
for a DPPC bilayer membrane and reported that 35-45 dyn/cm per
interface is needed to produce stable systems with reasonable values of surface area per molecule. Our result (35 dyn/cm per interface) agrees
well with their simulation.
Chain packing
Fig. 11 shows two-dimensional
radial distribution functions (RDFs) of the mass centers of DMPC
molecules (a), alkyl chains (b), and headgroups
(a). The positions of the first peak (~4.9 Å) observed in
the RDFs for the alkyl chains (Fig. 11 b) are consistent with the band observed in small-angle x-ray diffraction studies in many
hydrated phospholipid systems (Luzzati, 1968; Janiak et al., 1976). On
the other hand, the peak for the whole lipid is not clear (Fig. 11
a). These results indicate that alkyl chains are packed in
more or less regular arrays whereas whole lipids are not. The RDF for
the headgroup has a small peak at 4.4 Å, indicating that some
headgroup pairs interact strongly with each other (Pasenkiewicz-Gierula
et al., 1999).
View larger version (24K):
[in this window]
[in a new window]
|
FIGURE 11
(a) Two-dimensional radial distribution
function (RDF) for the mass center of the whole lipid (solid
line) and that of the head group (dotted line) projected on the
membrane plane. (b) Two-dimensional RDF for the mass
center of the alkyl chain (thick solid line).
Contributions from intra-lipid pairs (thin solid line)
and inter-lipid pairs (dotted line) are also shown
separately.
|
|
RDFs for the alkyl chains were calculated separately for those in the
same lipid and those in different lipids (Fig. 11 b). Such
information can be obtained only using a simulated membrane, and not
(easily) by experiments. The first peak at ~4.9 Å is clear in both
intralipid and interlipid RDFs, indicating that good packing occurs in
the spatial range of 4.9 Å regardless of whether alkyl chains belong
to the same lipid or different lipids. This reflects the importance of
alkyl chain packing for the integrity of the membrane.
Orientational correlation of shoulder vectors between neighboring
lipids
As described previously (see Fig. 6 and the related text), after
equilibration, no collective orientation of the DMPC shoulder vector
was found over 28 DMPC molecules (in one layer in system 1), and the
angular distribution of DMPC shoulder vectors about the molecular long
axis was rather uniform. However, this does not mean that there is no
short-range angular correlation of shoulder vectors (as they are
projected on the membrane plane) between neighboring DMPC molecules.
Because a DMPC molecule contains two alkyl chains, the overall shape of
the cross section of the hydrophobic region of the molecule is
anisotropic and is thought to be rectangular rather than circular. Such
anisotropy in the cross section of the hydrophobic domain could induce
orientational correlation between shoulder vectors of neighboring DMPC
molecules. This was examined, and the result is shown in Fig.
12.
View larger version (22K):
[in this window]
[in a new window]
|
FIGURE 12
Distributions of the angle made between neighboring
lipids' shoulder vectors. The vector that connects O21C3 projected
on the membrane surface is defined as the shoulder vector, which shows
the orientation of DMPC molecules about their long axes. Thick solid
line, first nearest neighbors; thin solid line, second nearest
neighbors; dotted line, third nearest neighbors.
|
|
Fig. 12 shows the distribution of the angle made between shoulder
vectors of neighboring lipids. Parallel and anti-parallel orientations
were not distinguished. The results in Fig. 12 show that the first
nearest neighbors prefer having parallel or anti-parallel orientations,
as expected from the shape of the cross section of the hydrophobic
region of the phospholipid. Such an orientational correlation extends
to the second nearest neighbors, but not to the third nearest
neighbors, showing that thermal effects overcome the geometrical effect
of each molecule in the distance to the second neighbors (~ 8 Å).
Although such an examination has not been possible in experimental
membranes, the development of a reliable computer model now enables us
to carry out such an investigation.
This result indicates that when the initial structures of computer
models are generated, merely randomizing the initial in-plane orientation of each molecule is not sufficient for equilibration in
terms of in-plane orientation. Sufficient time should be given for the
system to equilibrate so that it may develop the correct level of
orientational correlation between neighboring lipids.
Alkyl chain tilt
The average tilt angle of alkyl chains is nearly 0 (9°).
However, at the level of individual molecules, many are tilted. The distribution of alkyl chain tilt is shown in Fig.
13 a, which indicates that
the highest probability of the tilt angle is ~20° with respect to
the membrane normal. Some alkyl chains show tilt angles greater than
90°. During the production run (2.9-3.9 ns), 1-4 of 224 alkyl chains showed tilt angles greater than 90°. Representative
conformations of molecules showing such tilt angles are displayed in
Fig. 14. These chains are bent in the
middle and their methyl terminals reach the membrane surface, facing
the water layer. Such a conformation is consistent with a model that
was previously proposed to explain fast collisions of alkyl chain
methyl terminals with the membrane-water interface (Merkle et al.,
1987).
View larger version (20K):
[in this window]
[in a new window]
|
FIGURE 13
(a) Distributions of the alkyl chain
tilt angle. Alkyl chain orientation t is defined as the unit
vector pointing from the midpoint of C1 and C2 toward that of C13 and
C14, and alkyl chain tilt is defined as the angle between t
and the membrane normal (polar angle of t). (b)
Distribution of the angle between the y axis of the
simulation box and the projection of t on the membrane plane (azimuth
angle of t). (c) Projections of termini of
t on the membrane plane when their origins are brought
together to the center of the circle. Crosses (+) and small solid
circles ( ) represent and chains, respectively. Average
direction is denoted by a solid square ( ).
|
|
View larger version (58K):
[in this window]
[in a new window]
|
FIGURE 14
A snapshot of three lipid molecules (thick lines)
whose alkyl chains are bent in the middle while their terminals reach
the membrane surface in system 2 at 3.8 ns. Other lipid molecules are
omitted for clarity.
|
|
The angular distribution of t about the membrane normal is
shown in Fig. 13 b, which gives the angle of projection of t on the membrane plane with respect to the y
axis of the simulation box (which is in the membrane plane). Fig. 13
b shows that the angular distribution about the membrane
normal is almost uniform.
Fig. 13 c shows projections of the termini of t
vectors for all of the alkyl chains on the membrane plane when their origins are brought together to the center of the circle. As an alkyl
chain is tilted more, the tip of t moves outer from the
center of the circle. The average tip position of t is shown
by a solid square. This figure clearly shows that although the average
tilt angle is nearly 0, individual alkyl chains show a wide range of
tilt angles.
Correlation between alkyl chain tilt of neighboring chains
Although the average tilt angle is nearly 0, and the angular
distribution of t vectors about the molecular long axis is
broad, there must be short-range orientational correlations between
neighboring t vectors. The correlation of alkyl chain
orientation can be defined as
![⟨P<SUB>2</SUB>[<B><UP>t</UP></B><SUB><UP>i</UP></SUB>×<B><UP>t</UP></B><SUB><UP>j</UP></SUB>]⟩<SUB><UP>ensemble</UP></SUB>=½(3⟨<UP>cos</UP><SUP>2</SUP>&thgr;<SUB><UP>ij</UP></SUB>⟩<SUB><UP>ensemble</UP></SUB>−1),](/content/vol79/issue6/fulltext/3118/img006.gif)
|
(4)
|
where P2 is the second-rank
Legendre polynomial, ti is the unit
orientation vector of the ith chain, and
ij is the angle made by
ti and tj. For
the head groups, ti was defined as a unit
vector going from the phosphorus atom toward the nitrogen atom.
Inter-alkyl-chain distances are defined as the distances between the
mass centers of alkyl chains projected on the membrane surface.
Fig. 15 a shows the
correlations of orientations between alkyl chains and between
headgroups as a function of the distance between them. Two alkyl
chains, whether they are in the same molecule or belong to different
DMPC molecules, exhibited a high correlation at ~4.8 Å, which is the
same as the peak position for the radial distribution function of alkyl
chains (Fig. 11). These results indicate that neighboring alkyl chains
tend to exhibit similar tilt angles. Two headgroups that are in the
distance less than the first peak observed in the radial distribution
function (~4.4 Å, Fig. 11) also exhibit a high correlation. This
indicates that headgroup pairs interact strongly with each other
(Pasenkiewicz-Gierula et al., 1999). Correlation of alkyl chain
orientation persists long because of the overall alignment of alkyl
chains in the membrane, whereas no long-range correlation was observed
for headgroup orientations.
View larger version (21K):
[in this window]
[in a new window]
|
FIGURE 15
(a) Correlation between alkyl chain
orientations (tilt angles, thick solid line) and that of
headgroup orientations (broken line), plotted as a
function of the distance between mass centers. Correlations for alkyl
chains in the same DMPC molecule (thin solid line) and
those in different DMPC molecules (dotted line) were also evaluated
separately. Orientation is defined as the direction of the t
vector (see Eq. 4 and the text). The angle between pairs of
t vectors was obtained and
P2(cos) was calculated.
P2(cos) is shown as a function of
the distance between the mass centers (projected on the membrane plane)
of alkyl chains. (b) Distribution of the angle made between
neighboring alkyl chains. Thick line, alkyl chains within 7 Å distance; thin line, alkyl chains in the distance between 7 and 12 Å.
(c) Same as b, but normalized by sin to
exhibit the correct probability.
|
|
This result again indicates that when initial structures of computer
models are generated, merely randomizing the alkyl chain orientations
in the model is not sufficient. Such correlation should be examined to
test whether sufficient time was given for the system to equilibrate so
that the system may develop a correct level of orientational
correlation between neighboring lipids.
Fig. 15 b shows the distribution of the angle made between
neighboring alkyl chains. Based on the radial distribution function of
alkyl chain shown in Fig. 11, two alkyl chains within the distance of 7 Å (including the first peak in the RDF, thick line), and those in the
distance between 7 and 12 Å (including the second peak, thin line)
were used for examination of orientational correlation. The results in
Fig. 15 b indicate that more than half of neighboring chains
have tilt angles within 25°. Fig. 15 c shows the same
distribution but normalized by sin for the correction of the steric
angle. It indicates that the probability of finding well-aligned
t in neighboring alkyl chains is very high.
Chain conformation
For comparison of the Smol values
obtained here by the MD simulation with those obtained experimentally
by 2H-NMR and ESR spin labeling, we have to pay
special attention to the time scales with each spectroscopic method
(McConnell, 1976). Although an experimental result is the average over
many molecules (and thus close to the ensemble average), it should be
clearly understood that, in using these spectroscopic methods to
evaluate Smol values, the
spectroscopically measured Smol values are not the same as ensemble average values of
Smol values because the characteristic
time scales of the spectroscopic methods provide low-frequency cut
filters in the evaluation of Smol
values. The time scales for 2H-NMR and ESR spin
labeling are 10 µs and 10 ns, respectively, which differ by a factor
of 1000. On the other hand, MD simulation provides ensemble average
Smol values, which are expected to be the same as those that would be obtained by a hypothetical experimental method that had characteristic time scales longer than the time scales
of all dynamic molecular processes that affect
Smol values in the membrane.
As the characteristic time scale for 2H-NMR is
thought to be sufficiently long compared with the time scales of most
of the molecular processes that affect the measurement of
Smol values, including conformational
changes in DMPC and wobbling of the alkyl chains, time-averaging in
2H-NMR spectroscopy is not likely to cut off many
dynamic processes in the membrane. Therefore,
Smol values obtained by MD simulation were compared with those obtained by 2H-NMR, as
shown in Fig. 16 a. Because
the 2H-NMR data used here are for the -chain
in DPPC membranes obtained by Seelig and Seelig (1974),
Smol values from the MD data were calculated for -chains. Although the experimental data are for DPPC
(whereas our simulation uses DMPC), we think this comparison is
meaningful. Hubbell and McConnell (1971) showed that alkyl chain order
depends on the distance from the glycerol group. Furthermore, Kusumi et
al. (1986) directly compared the alkyl chain orders at different depths
in DMPC, DPPC, and DSPC membranes and found that the alkyl chain order
has the same dependence on the distance (carbon numbers) from the
glycerol group, irrespective of the alkyl chain length of PC that makes
up the membrane. The order parameter is calculated using the angle
between the instantaneous orientation of the segment vector,
ci, and the membrane normal (which is the
ensemble average orientation of ci). The
2H-NMR order parameter decreases slowly from the
segment C2 to C9. In this region, Smol
values of the MD-simulated membrane agree well with
Smol values by
2H-NMR. (The order parameters from the MD
simulation tend to be slightly less, suggesting that some
conformational changes or wobbling motion may occur at a rate slower
than 105
s1, which was the cutoff
in 2H-NMR evaluation of
Smol values. Although spectroscopic
techniques deal with great numbers of molecules, they cannot be
sensitive to processes that occur within time scales slower than the
spectroscopic time scales.)
and
TOP
ABSTRACT
INTRODUCTION
