INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS REFERENCES

The last decade has seen steady progress in the computer simulation of membrane-associated peptides and proteins (Merz and Roux, 1996; La Rocca et al., 1999; Forrest and Sansom, 2000). Given the experimental difficulties in studying these systems, it is certain that simulations will play an increasingly larger role in the elucidation of membrane structure and function.

Before simulation becomes an equal partner with experiment, however, there are numerous technical and theoretical problems that must be solved. One of these involves the limitations of periodic boundary conditions (PBC) when the relative number of lipids in the two leaflets of a bilayer requires adjustment during a simulation. In the usual implementation of PBC, denoted P1, unit cells are replicated by translation only. This is illustrated in the top panel of Fig. 1, which shows a particle in the primary cell (black) and its image (stippled) in P1 symmetry. If the black particle were to exit the primary cell from the right side of the top leaflet, the stippled particle would enter the cell at the top left. Edge effects are eliminated with this technique, and a microscopic patch of as few as 50-100 lipids can plausibly model the short time and length scale dynamics of a much larger bilayer.

View larger version (37K): [in this window] [in a new window]   FIGURE 1   Representation of a lipid bilayer in the three different periodic boundary conditions (PBC) discussed in this paper: P1 (top), P21 (middle), and Pc (bottom). The spheres represent lipid headgroups, with those in the top leaflet mostly white and one black, and in the bottom leaflet gray. The image of the black particle is stippled, and the arrows show the direction of exit and entry for these particles for each PBC.

Now consider the area/lipid (A_L) and number of lipids (N_L) in a microscopic patch of a membrane during the following hypothetical sequence: 1) insertion of a peptide into the "top" leaflet of the bilayer, parallel to the surface and with no disruption to the bottom leaflet; and 2) rearrangement of the peptide perpendicular to the surface, spanning the entire bilayer. In step 1, N_L is reduced and A_L is potentially changed in the top leaflet, but they are unaltered in the bottom leaflet. In step 2, N_L and A_L are the same in both leaflets, though different from the peptide free bilayer. In a cell membrane, lipids would flow into and out of the patch to relieve the local stress as the peptide inserts and rotates. Simulation of this sequence using P1 PBC would require an explicit, and largely ad hoc, adjustment of lipids at each step. Even the simpler case of inserting a peptide into a bilayer for simulation is typically carried out using approximate methods relying on estimates of projected surface areas (Woolf and Roux, 1994; Berneche et al., 1988; Tang et al., 1999). In principle, it is possible to insert or delete lipids during molecular dynamics by simulating in the grand canonical ensemble (Lynch and Pettitt, 1997), though such an approach is not currently practical for membranes. Realistic simulation of peptide insertion and rearrangement, and evaluation of the free energy changes involved in these processes, is essential for an understanding of the interaction of antimicrobial and lytic peptides with membranes (Shai, 1999; Bechinger, 1999).

This paper considers two alternative transformations in PBC designed to allow exchange of lipids between leaflets. These are also sketched in Fig. 1. In each case, a particle exiting from the top leaflet of the primary cell enters the bottom leaflet. Fig. 2 illustrates the transformations involved. In P2₁, more completely denoted P2₁(x, y + 1/2, z), the primary cell is rotated about the y axis and translated along the y axis by one-half a cell length; the rotation has the effect of transforming x to x and z to z. In a simulation of a lipid bilayer, the effect of the P2₁ boundary condition is that a lipid exiting either leaflet re-enters the opposite leaflet via an orthogonal face. In Pc(x + 1/2, y + 1/2, z), the cell is inverted through a mirror plane at z = 0. P1(x, y, z), by contrast, is translationally invariant. No lipid has a close interaction with its symmetry-related image in any of the three preceding boundary conditions.

View larger version (35K): [in this window] [in a new window]   FIGURE 2   The P21(x, y + 1/2, z) and Pc(x + 1/2, y +1/2, z) transformations, looking down the normal (z axis) to the bilayer (xy plane). The primary unit cell is white, and located in the center; the four particles in the corners of the cell reside on the top leaflet. In the adjacent cells (dark gray), the same four particles are in the bottom leaflet, are rotated with respect to the primary cell for P21, and inverted for Pc. The outer light gray cells are the result of a second transformation, and the four particles return to the top leaflet.

Molecular dynamics (MD) simulations of five different systems are reported. These include 1) bulk water; 2) a water/vapor interface; 3) a monolayer film of octane initially localized to a single face on a slab of water; 4) bilayers of dipalmitoyl phosphatidylcholine (DPPC) with equal numbers of lipids per side, and with an initial mismatch consisting of 40 lipids on one leaflet and 32 on the other; and 5) a dioleoyl phosphatidylcholine (DOPC) bilayer containing a monomer of melittin in one leaflet. All systems except (3) were simulated in both P1 and P2₁, and the first two were also simulated with Pc. As appropriate and as described in the Methods section, ensembles included constant volume, constant surface area, and constant surface tension. In addition to thoroughly testing the new method and illustrating the time scale of rearrangements in bilayers, these simulations highlight the limitations in studying membranes with only P1 boundary conditions.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS REFERENCES

General considerations

All MD simulations were carried out with the program CHARMM (Brooks et al., 1983) with the CHARMM 27 academic all-atom parameter set (Feller and MacKerell, 2000) and a modified TIP3P water (Jorgensen et al., 1983; Durell et al., 1994). Electrostatic interactions were calculated via the particle mesh Ewald (PME) summation (Essmann et al., 1995) using a real space cutoff distance of 10 Å (except as noted below), a screening parameter  = 0.34 Å¹, an fft grid density of ~1 Ź, and a fifth-order polynomial spline for interpolation between the grid points. For the bulk water, octane film, and matched bilayers, the Lennard-Jones (LJ) 6-12 potential was switched smoothly to zero over the region from 8 to 10 Å. For the mismatched DPPC bilayers and DOPC/melittin systems the LJ potential was switched to zero over the region from 7 to 11 Å, and the real space cutoff for Ewald summation was set to 11 Å. Nonbond and image lists were updated every 20 integration steps. Net translation and rotation of the system were removed every 1000 steps. All bonds involving hydrogen were fixed at their equilibration distances using the SHAKE algorithm (Ryckaert et al., 1977). A time step of 1 fs was taken with a modified leap-frog Verlet integration technique.

All three PBC, P1(x, y, z), P2₁(x, y + 1/2, z), and Pc(x + 1/2, y + 1/2, z) were generated using the crystal facility available in CHARMM; i.e., existing code was sufficient. Simulation cells were either cubic (cell lengths a = b = c, and angles  =  =  = 90°) or tetragonal (a = b  c,  =  =  = 90°); c will be taken to be coincident with the z axis, and normal to the interface for the water/vapor, octane film, and bilayer systems. Systems initially simulated in P1 were centered at the origin with a and b coincident with x and y. For the alternate symmetries with the same number of atoms in the asymmetric unit, the cells were rotated by 45° in the xy plane and translated by a/4 in the +x direction (see Fig. 2). The center of the bilayer or water slab was positioned at z = 0; care must be taken to ensure this, or a molecule will be badly out of register upon rotation about the y axis or z axis inversion. With these placements, P2₁ and Pc transformations lead to the desired up/down checkerboard arrangement of adjacent cells (Fig. 2). Note that the 45° rotation by itself is insufficient; the cells would overlap after the transformation without the x axis offset. As a practical matter, atomic coordinates were scaled by 0.98 when converting between symmetries, while maintaining the unit cell dimensions. The systems were then minimized to establish the new contact surfaces and restore appropriate interatomic spacing before further simulation.

Simulations using P2₁ or Pc symmetries and PME required 20-30% more cpu time than those carried out in P1, depending on system size. This is because the PME method as currently implemented computes the grid based on a full unit cell, thus doubling the number of grid points for the two asymmetric unit systems. Note also that the chirality of a particle is reversed when changing leaflets in Pc. Consequently, this symmetry should only be used for nonchiral molecules or for systems representing a 50:50 chiral mixture.

Simulations were carried out in a variety of ensembles to test the different boundary conditions thoroughly. These included NVE, NVT, NPAT, and NPT, where N is particle number, V is volume, E is energy, T is temperature, P is normal pressure, A is surface area, and  is surface tension. The simulation of interfaces and membranes using these ensembles has already been described (Zhang et al., 1995; Feller et al., 1995a); numerical values for parameters are listed below for each system. Unless it is sufficiently clear, a simulation will be denoted by both the ensemble and boundary condition, e.g., NPAT/P2₁.

Finally, we consider the question, are there one or two surfaces? Formally, there is only a single surface, like a Mobius strip. Whether this simulated asymmetric unit is represented as a bilayer or as a larger single surface depends on the image centering procedure that is used, but this choice will not affect the trajectory. For ease of exposition we consider a bilayer of 36 lipids/leaflet. Consider Fig. 2 and suppose that there are 36 lipids on average in the top leaflet (light squares) and 36 in the lower leaflet (gray squares). The entire bilayer could be represented as a single rectangular surface consisting of a light square and dark square with 72 lipids. A lipid may diffuse to every position of either surface without ever passing through the interior. There will be density fluctuations which, in effect, are manifested in particle exchange between leaflets (e.g., 38 lipids in the upper and 34 in the lower).

There is a difficulty, however, when counting the number of lipids on each side. For example, for purposes of placement in a periodic system the symmetry location of a molecule might be specified by its center of mass. The symmetry transformation that results in the shortest distance to a specified point (usually the origin, though not for the P2₁ case) is chosen on periodic intervals. In CHARMM, the location is determined by the center of volume of a rectangular prism enclosing all of the atoms of the molecule. Either image-centering method leads to integer counts. However, a lipid near an edge can have atoms (image or primary) in both leaflets, and different image-centering procedures can lead to different placements. Therefore, it is potentially misleading to calculate quantities such as surface area per lipid per leaflet for systems in P2₁ or Pc symmetry using simple integer counts, though trends may be ascertained.

Bulk water and water/vapor interface

A previously equilibrated box of 1340 water molecules with dimensions of 34.2 Å × 34.2 Å × 34.2 Å was used to initially set up the new PBC. This bulk system was then used to construct a water-vacuum interface by increasing the length of the box in the z direction from 34.2 Å to 64 Å, creating a vacuum volume with a distance of ~30 Å between the water/vacuum planar interfaces.

NVE and NVT MD simulations were carried out on both systems under P1, P2₁, and Pc to establish feasibility. Simulations were 210 ps in length, with the first 10 ps used for equilibration. For NVT simulations, the temperature was maintained at 293 K using a piston (Hoover, 1985) with an effective mass of 20,000 kcal/ps². Simulations were carried out using two, four, or eight processors of the National Institutes of Health LoBoS Linux cluster (http://www.lobos.nih.gov).

Octane film

The initial conditions for this system consisted of 24 octanes on one side of a slab of 1340 waters (see Fig. 4, left). Ten statistically independent starting structures were obtained by extracting coordinates from a previous simulation trajectory of a 30 Å slab of octane between two water layers and then deleting all octanes beyond ~6.5 Å from one of the water/octane interfaces. (Extensive discussion of simulation of the water/octane system is included in Zhang et al. (1995) and Feller et al. (1996).) NVT/P2₁ simulations of 210 ps each were run on these 10 independent starting structures; one of these 10 was also simulated in NVT/Pc. The cell dimensions were identical to those of the preceding water/vacuum system.

DPPC bilayers

Initial test systems were a bilayer with 72 DPPC and 2094 waters, and a bilayer with 36 DPPC and 1047 waters. The former was extracted from the final (800 ps) point of a previous DPPC simulation (Feller et al., 1997a). The latter was constructed by dividing the 72-lipid system in half. The 72- and 36-lipid systems were simulated for 200 and 400 ps at 323 K, respectively, in both NPAT/P1 and NPAT/P2₁.

DPPC mismatch

Twelve coordinate sets consisting of 36 lipids/leaflet were randomly picked from a previous 10 ns DPPC bilayer simulation (Feller et al., 1999). Four lipids for each set were then rotated 180° about the Cartesian x or y axes to yield bilayers with 40 lipids in one leaflet and 32 in the other (both the lipids and the axis of rotation were randomly chosen). The systems were stabilized by minimization using fixed constraints on the lipid just rotated, and harmonic restraints in z for water hydrating the leaflet that just lost a lipid. A second stage of minimization was performed with the fixed constraints removed, and a third stage with no harmonic restraints. Each minimization stage comprised 50 steps of the steepest descent (SD) algorithm, followed by 200 steps of the adopted basis Newton-Raphson (ABNR) algorithm. An MD simulation with heating from 223 K to 323 K over 1 ps at NPAT/P1 completed the stabilization of each rotated lipid. After examination of energy, pressure, and cell dimensions, six were randomly chosen for further study.

Each of the six models was first converted from P1 to P2₁ symmetry by applying the z axis rotation, x axis translation, and coordinate scaling noted above. Each was then minimized for 50 steps of the SD algorithm, followed by 500 steps of the ABNR algorithm, with image centering enabled only for the water molecules, to prevent any lipid leaflet exchange during the minimization, and then simulated for 400 ps of MD in NPAT/P2₁.

High-hydration DOPC

A fully hydrated DOPC bilayer (2358 waters for 72 lipids, or 32.75 waters/lipid) at 303 K was constructed in stages by first adding additional water to a low-hydration DOPC coordinate set (Feller et al., 1997b) and equilibrating the system with 400 ps of MD under NPAT/P1 conditions. The ensemble was then switched to NPT to expand the surface area from the low hydration value of 59.3 Å²/lipid to the experimental value of 72.5 Å²/lipid at a hydration of 32.8 waters/lipid (Nagle and Tristram-Nagle, 2000). The surface tension was first set to 50 dyn/cm, changed to 75 dyn/cm at 100 ps to accelerate the expansion, and then back to 50 dyn/cm at 200 ps to decelerate as the experimental target area was approached and finally attained at 208 ps. The system was further equilibrated for 400 ps at NPAT. The final coordinate set of this simulation was converted to P2₁ via a minimization protocol of 50 and 200 steps of SD and ABNR algorithms, respectively. An additional 400 ps NPAT/P2₁ simulation was carried out, and the final coordinate set used to create the DOPC/melittin system described below.

DOPC/melittin

Several models of the DOPC/melittin complex with varying numbers of lipids removed and different methods for insertion were investigated, and the comparisons will be presented elsewhere. Results for a single illustrative case starting with the melittin structure obtained from the neutron diffraction studies of White and co-workers (Hristova et al., 2001) will be described here.

Melittin (coordinate set B = 152 from Hristova et al., 2001) was inserted by model building into one leaflet of the DOPC bilayer parallel to the interface plane at the approximate level of the glycerol and with the hydrophobic face pointed toward the bilayer center, as indicated by experiment. Six lipids and 50 overlapping waters were then deleted (leaving 30 lipids in the melittin-containing leaflet, 36 in the other, and 2308 waters), and 6 chloride ions were added to attain electroneutrality. This lipid distribution was comparable to the initial condition constructed by Berneche et al. (1998) of a melittin in a dimyristoyl phosphatidylcholine (DMPC) bilayer: 17 lipids surrounding melittin and 24 in the lower leaflet. Given that the surface areas per lipid of pure DMPC and DOPC at full hydration are 59.6 and 72.5 Å², respectively (Nagle and Tristram-Nagle, 2000), the area of 7 DMPC is close to that of 6 DOPC.

The following minimization protocol was then carried out in P2₁ at fixed surface area: 1) 20-step SD algorithm with LJ term only and fixed constraints were placed on peptide atoms; 2) 20-step SD algorithm, 50-step ABNR algorithm with full potential restored; 3) 20-step SD algorithm, 50-step ABNR algorithm with the fixed constraints on peptide atoms removed and harmonic restraints placed on the C positions; 4) 20-step SD algorithm, 50-step ABNR algorithm with all restraints removed. The system was then simulated for 450 ps at NPT/P2₁, with  = 25 dyn/cm and helix backbone harmonic shape restraints applied for the first 200 ps.

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS REFERENCES

Bulk water and water/vacuum interface

A variety of equilibrium properties were evaluated for bulk water and the water/vacuum interface for the P1, P2₁, and Pc periodic boundary conditions. Radial distribution functions for the bulk (data not shown) and density profiles for the interface (Fig. 3) are close to superimposable for the three boundary conditions, and there are no statistically significant differences in the surface tension, isotropic pressure, or energy (Table 1). Surface tensions were statistically equivalent to zero for the bulk, as desired, and to 53 dyn/cm for the interface, as calculated earlier for the TIP3 model (Feller et al., 1996); the slightly negative isotropic pressures for the bulk are characteristic of the model. The total Hamiltonians are also listed in Table 1. Their low standard error indicates that the fluctuations were small and the trajectories were stable (the averages are not statistically equivalent because the systems were not prepared identically). Because results of these preliminary simulations were positive, more complex systems were investigated.

View larger version (18K): [in this window] [in a new window]   FIGURE 3   Density distribution function normal to the interface, (z), from NVT simulations of the water/vacuum interface under P1 (solid line), P21 (dashed line), and Pc (dotted line), periodic boundary conditions.

View this table: [in this window] [in a new window]   TABLE 1   NVT simulation averages and standard errors for surface tension , isotropic pressure P, total energy E, and total Hamiltonian for bulk water and water/vacuum under three different periodic boundary conditions (PBC)

Octane film

As shown in Fig. 4 for one of the 10 systems run in P2₁, individual octanes readily diffuse to the opposite face of the water slab from their initial localized placement. On average, the present system took ~100 ps to equilibrate (Fig. 5, thick line), though fluctuations were substantial (Fig. 5, thin lines). Averaging over the second 100 ps, the number of octanes on the top face are 5.6, 7.0, 9.8, 11.5, 12.7, 13.8, 16.0, 16.7, 17.9, 19.1. These (sorted) data yield an average of 13.0, a standard deviation of 4.6 and standard error of 1.4. Hence, the average is statistically equivalent to 12 per side. It is useful to compare these results with those expected from a binomial distribution with parameters n (the number of Bernoulli trials, or octanes) = 24 and p (the probability of an octane being in the top surface) = 0.5. The standard deviation of such a distribution (DeGroot, 1975) is p  2.4, which is significantly smaller than the observed standard deviation of 4.6. This result implies that the individual octanes are not independent. Rather, they tend to aggregate in a lens on the water surface.

View larger version (67K): [in this window] [in a new window]   FIGURE 4   Three snapshots from a simulation of an octane monolayer on a slab of water using P21 periodic boundary conditions: t = 0 ps (left), 20 ps (middle), and 200 ps (right).

View larger version (41K): [in this window] [in a new window]   FIGURE 5   Average number of octane molecules in the top half of the water slab versus time (thick line). Instantaneous counts for each of the 10 trajectories are also included (thin lines). All simulations began with 24 octanes in the top half.

A single simulation was also carried out with Pc, with qualitatively similar results. Because the chirality reversal inherent in this PBC precludes its application to most biophysical systems, the remaining simulations presented here were carried out in either P1 or P2₁.

Lastly, while equilibration between sides of the water slab took place almost exclusively along the surface, occasionally an octane crossed to the other side via the vacuum space. Hence, slow mixing would have been observed had the simulation been carried out with P1.

Pure DPPC bilayers

Simulations of 36- and 72-lipid systems were stable in both P1 and P2₁, and calculated surface tensions for all four systems were statistically equivalent. This equivalence would not necessarily hold for all bilayer sizes. Undulatory modes contribute to bilayer surface tensions (Feller and Pastor, 1996, 1999), and the modes of the P1 and P2₁ systems are not identical. The present results, however, are consistent with our earlier work (Feller and Pastor, 1996), where bilayers of 16 and 36 lipids per leaflet had similar surface tensions when simulated in P1. Hence, collective modes make a relatively small contribution to surface tension in bilayers of these sizes, and the P1 and P2₁ systems have similar surface tensions. The deuterium order parameters were also evaluated. Table 2 lists representative values, including the average for carbons 4-8, a common measure of the order of the plateau region (Feller et al., 1997a). The differences in values for each system are well within their standard errors, indicating equivalence of the methods.

View this table: [in this window] [in a new window]   TABLE 2   Surface tensions and selected deuterium order parameters (SCD) with standard errors for the 4 DPPC (matched) bilayer simulations. C4-C8 denotes the average for carbons 4-8

Fig. 6 plots the equilibration of lipids during the mismatch simulations, where the initial configurations consisted of 40 lipids in the top leaflet and 32 in the bottom. There is a rapid drop to 37 lipids in the top within 10 ps, averaged over the six trajectories. This is followed by sometimes substantial fluctuations (as high as 41 and as low as 33), and the average is still close to 37 (rather than 36) after 400 ps.

View larger version (36K): [in this window] [in a new window]   FIGURE 6   Average number of lipids in the top leaflet of the DPPC bilayer versus time (thick line), and instantaneous counts for each of three of the six trajectories (assorted thin lines). All simulations began with 40 lipids in the top leaflet and 32 in the bottom.

The preceding bilayer simulations were carried out at constant particle number, normal pressure, surface area, and temperature (NPAT) ensemble. Consequently, the natural mechanism for reducing stress in the mismatched systems is lipid exchange between the leaflets. The mismatched systems were also simulated in the constant surface tension ensemble, NPT, with  = 20 dyn/cm, and with damping on the pressure piston (Feller et al., 1995b). These simulations showed a rapid expansion in the surface area, followed by a somewhat slower exchange of lipids and a contraction of surface area.

Melittin in DOPC

Melittin was inserted into the fully hydrated and expanded bilayer as described in the Methods. The central result of this section is shown in Fig. 7: within 200 ps the number of lipids in the bottom (melittin-free) leaflet decreases from an initial value of 36 to 32 ± 1. This was accompanied by a small contraction of 2.6% in the surface area. The observed flow of lipids into the melittin-containing leaflet highlights how potential inaccuracies in initial model building can be overcome by simulation in P2₁. Fig. 8 shows the final 450 ps coordinate frame. Averages for melittin calculated over the final 200 ps show relatively little change from the experimental structure model with B = 152: a root-mean-square deviation (rmsd) of 1.98 Å; average z position of 16.8 ± 0.1 Å compared with 17.1 ± 0.3 Å; and radius of gyration (R_g) of 13.04 ± 0.03 Å compared with 13.55 Å. Difference in the simulated values and the B = 152 model are comparable to those between that model and another model with B = 175 derived by Hristova et al. (2001) from their diffraction data. The simulated values do differ significantly (rmsd = 3.82) from the x-ray structure of melittin where R_g = 11.7 Å (Terwillinger and Eisenberg, 1982). Lastly, as evidenced by a comparison of calculated order parameters from bilayers of melittin/DOPC and pure DOPC (Fig. 9), melittin significantly perturbs the lipids near the top of the chains, yet appears to leave the lower carbons relatively unchanged. This behavior is qualitatively similar to that observed experimentally in DPPC/melittin bilayers by Dufourc et al. (1986).

View larger version (24K): [in this window] [in a new window]   FIGURE 7   The number of lipids in the bottom leaflet of a DOPC/melittin bilayer versus time. The simulation began with 30 lipids and a single melittin in the top leaflet, and 36 lipids in the bottom.

View larger version (69K): [in this window] [in a new window]   FIGURE 8   Top down (top) and side (bottom) views of melittin in a DOPC bilayer from the final (450 ps) frame of the simulation in P21 symmetry. Melittin in space-filling representation with heavy atoms of hydrophobic residues are yellow and other carbons gray, nitrogens blue, and oxygens red; lipid carbons are gray, oxygens red, phosphate groups green, and the quaternary amine group purple; waters are blue and chloride ions large green.

View larger version (17K): [in this window] [in a new window]   FIGURE 9   Calculated deuterium order parameters of lipid chains from simulation of pure DOPC (squares) and DOPC/melittin (circles).

DISCUSSION AND CONCLUSIONS

We have demonstrated that the P2₁ and Pc boundary conditions yield equivalent results to P1 for bulk water and a water/vacuum interface (Fig. 3 and Table 1), and allow efficient spreading of octane film (Figs. 4 and 5). P2₁ was further tested on pure DPPC bilayers with an equal number of lipids per leaflet (Table 2), and ones with an initial mismatch in lipid number (Fig. 6). Lastly, a simulation of a DOPC/melittin membrane in P2₁ showed equilibration from a stressed initial condition (Fig. 7).

These are not the first MD studies of fluid systems with nonstandard boundary conditions. Hyvonen et al. (1997) simulated a bilayer with S₄ symmetry (the leaflet was flipped 180° and rotated by 90° about the axis normal to the interface). Wong and Pettitt (2000) simulated a silicate/water interface in the Pb spacegroup (the lower half of the unit cell is obtained by a reflection and translation). In both of these studies the unit cell contains only a single interface, and the methods were proposed as a means to reduce computer cost relative to a full P1 simulation with twice the number of atoms. The use of P2₁ presented here results in the same cost savings. However, we justify its use as a means to equalize the chemical potential on both leaflets of the bilayer; i.e., redistribute lipids between leaflets. In contrast to the P2₁, the S₄ symmetry has symmetric image atoms proximal to the corresponding primary atoms, resulting in undesirable correlation between leaflets. This correlation is especially problematic when simulating imbedded proteins. As in the case of Pc, the Pb symmetry reverses chirality, making it inappropriate for most biologically interesting membranes. Other space groups beyond those discussed here are also possible candidates, though the P2₁ appears to be optimal for bilayers.

Because simulations carried out in P2₁ and Pc can require somewhat more computer time than those at P1 depending on the cutoff method, it may be reasonable in some applications to equilibrate in the former and carry out the "production simulation" in the latter. There are two ways of doing this. One is to maintain the number of atoms; i.e., simulate in the same sized P1 system. This will require a period of re-equilibration to relax bad contacts between newly formed periodic edges. The other is to double the number of atoms by generating the other asymmetric unit and including it in a larger P1 system. This will not require an intermediate equilibration. It is also advisable to simulate in a constant pressure or constant surface tension ensemble, at least during equilibration, and to apply stochastic damping to the pressure pistons at this stage. Damping can be removed at the later stages of the simulation. Differences between P1 and P2₁ in simulations of lipid bilayers, possible in principle because coupling between the two leaflets is stronger in P2₁, do not appear to be significant for the systems studied here (18 and 36 lipids/leaflet).

As shown in DPPC bilayers, there appear to be at least two time scales of lipid equilibration in response to a mismatch: a rapid redistribution within ~10 ps and a much slower decay, at least hundreds of ps (Fig. 6). Nevertheless, an ensemble of reasonably equilibrated states was reached from the clearly stressed initial conditions, and the time scale is considerably shorter that that of lateral diffusion of the lipids (Pastor and Feller, 1996). As is consistent with small size systems where fluctuations in the surface area per lipid are large (Feller and Pastor, 1999), fluctuations in the number of lipids in a leaflet are also large. These natural fluctuations in the number of lipids in a small patch of bilayer could be related to formation of pores and other defects.

A DOPC/melittin system was constructed with 30 lipids in the melittin-containing (top) leaflet and 36 in the other. Importantly, subsequent simulation in P2₁ showed a net migration of four lipids into the melittin-containing leaflet (Fig. 7). It therefore appears that significantly fewer than six lipids needed to be removed because melittin condenses DOPC. The present result, perhaps related to the membrane lytic function of melittin, should be regarded as preliminary. Note that the final number of lipids in the lower leaflet is relatively small, only 32; i.e., there are more lipids in the melittin-containing leaflet. It is not yet clear, for example, whether lipids would continue to condense around the peptide or whether the bilayer would expand if the number in the lower leaflet is increased to 36 (approximately the number expected for a pure DOPC bilayer). Hence, it is not only important to allow redistribution of lipids between leaflets, but one must be able to compare the free energies of different configurations of the different alternatives.

The new boundary conditions do not substantially change the underlying dynamics, so one expects that systems close to equilibration would relax in a similar manner in all boundary conditions discussed here. The main advantage of P2₁ is in systems where there is no route to a satisfactory equilibrium structure in P1; i.e., the lipids must switch sides, as in mismatched DPPC and melittin/DOPC bilayers. Other systems that could be simulated to advantage in P2₁ include membrane proteins with noncylindrical shapes and bilayers with high curvature, such as vesicles or fusion intermediates. In closing, we anticipate that the P2₁ boundary condition will be exceedingly useful for simulations of complex membranes.