21

Membrane Simulations

Douglas J. Tobias

University of California at Irvine, Irvine, California

I.INTRODUCTION

Biological membranes hold cells together and divide them into compartments. They are also home to the lipid-soluble proteins that comprise roughly one-third of the genome and perform a variety of functions, including energy production and storage, signal transduction, and the formation of channels for the transport of substances into and out of cells and between their compartments. In 1972, Singer and Nicolson [1] assimilated available data into a model for the molecular organization of biological membranes that is now taken for granted. According to this ‘‘fluid mosaic’’ model, a lipid bilayer forms the membrane matrix in which proteins are embedded. This matrix is fluid at physiological temperatures, and both lipids and proteins are free to diffuse in the plane of the membrane.

Biological membranes contain a complex mixture of several different types of lipids. Membrane lipids are amphiphilic molecules, with both polar and nonpolar substituents. A common membrane lipid that has been the subject of numerous experimental and theoretical studies, dipalmitoylphosphatidylcholine (DPPC), is diagrammed in Figure 1. Like other membrane lipids containing two chains, DPPC is classified as an insoluble, swelling amphiphile [2]. These lipids exhibit a variety of liquid crystalline phases that may be lamellar, cubic, or hexagonal, depending on the temperature and water content. Lamellar phases consist of stacks of lipid bilayers separated by layers of water (Fig. 2). The structures and mechanical properties of lipid bilayers depend on the amount of water present. As water is added to a lamellar phase, the spacing between the bilayers increases until the ‘‘fully hydrated’’ state is reached. At this point, which in DPPC is about 28 water molecules per lipid, additional water molecules do not go between the bilayers but rather go into a bulk water phase that is distinct from the lamellar bilayer water phase.

Hydrated bilayers containing one or more lipid components are commonly employed as models for biological membranes. These model systems exhibit a multiplicity of structural ‘‘phases’’ that are not observed in biological membranes. In the state that is analogous to fluid biological membranes, the ‘‘liquid crystal’’ or Lα bilayer phase present above the main bilayer phase transition temperature, Tm, the lipid hydrocarbon chains are conformationally disordered and fluid (‘‘melted’’), and the lipids diffuse in the plane of the bilayer. At temperatures well below Tm, hydrated bilayers exist in the ‘‘gel,’’ or Lβ, state in which the mostly all-trans chains are collectively tilted and pack in a regular two-dimensional

465

466

Tobias

Figure 1 Chemical structure and space-filling representation of a phosphatidylcholine, DPPC. Different parts of the molecule are referred to by the labels at the left; together the choline and phosphate are referred to as the headgroup, which is zwitterionic. In the space-filling model, H atoms are white, O and P gray, and C black. (From Ref. 55.)

array. Because the biologically relevant liquid crystalline phase is highly disordered, it is not amenable to atomic resolution structure determination. However, electron density and scattering length density profiles as a function of depth in the membrane have been mapped out by X-ray and neutron diffraction [3–6], and deuterium NMR has been used to characterize the average order in the hydrocarbon chains [7]. The dynamics of lipids and water molecules in membranes on time scales presently accessible to simulations have been characterized by various experimental techniques, including neutron scattering [8,9] and NMR [8,10].

In recent years molecular dynamics (MD) simulations have emerged as a useful tool for filling in some of the details on the molecular structure and dynamics of membranes that are not available from experiments and for providing feedback on models used to interpret experimental data. Computer hardware and simulation methodology have matured to the point where stable simulations of membranes can now be routinely performed for at least a few nanoseconds, and many of the results of membrane simulations agree remarkably well with available experimental data [11]. With the help of simulations, Singer and Nicolson’s cartoon of a membrane drawn as a flat bilayer of lollypop-like lipids is evolving into images such as the one shown in Figure 2.

This chapter describes some of the technical aspects of simulating membranes, presents results that illustrate the novel insight into membrane structure and dynamics that can be provided by simulation, and discusses the correspondence of the emerging atomic scale picture with the results of NMR and X-ray and neutron scattering experiments. We restrict our attention to pure lipid bilayers, which have largely been the focus of the field of membrane simulations to date (recent reviews can be found in Refs. 11 and 12), with an emphasis on dynamics, which has been somewhat neglected compared to structure. However, we should mention that although the machinery for simulating membranes has been developed and tested on pure lipid bilayers, the number of applications to more

Membrane Simulations

467

Figure 2 Snapshot from an MD simulation of a multilamellar liquid crystalline phase DPPC bilayer. Water molecules are colored white, lipid polar groups gray, and lipid hydrocarbon chains black. The central simulation cell containing 64 DPPC and 1792 water molecules, outlined in the upper left portion of the figure, is shown along with seven replicas generated by the periodic boundary conditions. (From Ref. 55.)

complicated systems has been growing rapidly. A by no means exhaustive list includes the transport of small molecules (water [13], benzene [14], anesthetics [15,16]) and ions [17] within and across bilayers, cholesterol effects on bilayers [18], and membrane-bound peptides and proteins (surface-bound tripeptides [19,20], transmembrane polyalanine [21], an amphipathic helical peptide [22], an amphipathic segment from the human corticotrophin releasing factor [23], the lytic peptide mellitin [24], the antimicrobial peptide dermaseptin B [25], bacteriophage Pf1 coat protein [26], the transmembrane domain of ErbB-2 tyrosine kinase receptor [27], and bacteriorhodopsin [28–30]), including pores and ion channels (gramicidin [31–34], E. coli OmpF porin [35], alamethicin [36], and the influenza A M2 channel [37]).

II. MOLECULAR DYNAMICS SIMULATIONS OF MEMBRANES

A.System Size and Construction

Algorithms and computer hardware have limited the size of most of the all-atom membrane simulations performed to date to 50–100 lipids plus water. However, with the increasing availability of parallel computers, larger systems containing up to 1000 lipids are starting to be simulated [75], and systems containing more than 100 lipids will soon be routine.

468

Tobias

In the typical setup, the lipids are arranged in a bilayer, with water molecules on both sides, in a central simulation cell, or ‘‘box,’’ which is then replicated by using threedimensional periodic boundary conditions to produce an infinite multilamellar system (Fig. 2). It is important to note that the size of the central cell places an upper bound on the wavelength of fluctuations that can be supported by the system.

Because there are no crystal structures of fully hydrated lipid bilayers to start from, the initial configurations for the first membrane simulations had to be constructed ‘‘by hand.’’ Many of the early simulations were based on the crystal structure of dimyristoylphosphatidylcholine dihydrate [38], in which the DMPC molecules form a bilayer with slightly tilted all-trans hydrocarbon chains. The unit cell was replicated to give the desired number of molecules, and the lipid center-of-mass positions were scaled to give the desired area per lipid. Then the system was hydrated and equilibrated to allow hydrocarbon chains to melt and the water molecules to penetrate through the headgroup regions. A less used approach that reduced the equilibration time was to pack together lipids randomly chosen from a library of configurations generated by simulating a single lipid in a mean field [39]. Woolf and Roux extended this approach by using prehydrated lipids [31]. Presently, most simulations are initiated from the end of a previous simulation, and lipid bilayer coordinates are generally provided upon request by authors of simulation papers. Incorporating solutes into bilayers can be tedious. The simplest approach, which consists of simply deleting enough lipids to create ample space, works well in some cases (e.g., for inserting cholesterol molecules or transmembrane helices), but more sophisticated approaches are more generally useful. The latter include using radial forces to create spherical or cylindrical holes in the bilayer [36], and growing groups of soft spheres at locations that create a cavity in the membrane with the desired size and shape [26].

B. Force Fields

The force fields used in atomistic membrane simulations have the same form as the analytical, empirical molecular mechanics potentials used in classical simulations of proteins and nucleic acids (see Chapter 2). These potential functions contain harmonic terms for deforming bonds and bond angles, periodic and harmonic potentials for torsions, and van der Waals and Coulomb nonbonded interactions. The atomic charges are often obtained from quantum chemical calculations on lipid fragments. Most of the other potential parameters are taken from existing force fields for proteins and nucleic acids, which are generally built up from parametrizations of model compounds. Details on biomolecular force field parametrization are given in Refs. 40 and 41, and in Chapter 2. The quality of lipid force fields may be evaluated with simulations of crystals of lipid fragments [40,42,43] and by checking their ability to reproduce well-established experimental results in simulations of hydrated lipid bilayers (see Section III.A). The level of accuracy that we have been able to achieve is typified by the densities of phospholipid and cholesterol crystals from constantpressure MD simulations with fully flexible unit cells, plotted in Figure 3 versus the experimental values.

Calculation of the energies and forces due to the long-range Coulomb interactions between charged atoms is a major problem in simulations of biological molecules (see Chapter 5). In an isolated system the number of these interactions is proportional to N 2, where N is the number of charged atoms, and the evaluation of the electrostatic interactions quickly becomes intractable as the system size is increased. Moreover, when periodic

Membrane Simulations

469

Figure 3 Comparison of the densities (in g/cm3) of model compounds for membrane lipids computed from constant-pressure MD simulations with the corresponding experimental values. The model compounds include solid octane and tricosane, liquid butane, octane, tetradecane, and eicosane, and the glycerylphosphorylcholine, cyclopentylphosphorylcholine monohydrate, dilaurolyglycerol, anhydrous cholesterol, cholesterol monohydrate, and cholesterol acetate crystals. (Models from Refs. 18, 42, and 43).

boundary conditions are used, the direct sum over all the Coulomb interactions in the periodic system is conditionally convergent. There are two common solutions to these problems. The first, which was used in most biomolecular simulations until recently, is to employ a spherical truncation (‘‘cutoff’’) scheme where the electrostatic interactions

˚

are smoothly switched off at a cutoff distance (typically around 10 A). When used in conjunction with a neighbor list, the computational work is proportional to N. An alternative approach, which is essentially exact for crystals, is to include all of the electrostatic interactions by using the Ewald method (or a related technique) for summing long-range interactions in a rapidly convergent fashion in an infinite periodic system. The disadvantages of the Ewald method are that it might enhance the artificial periodicity of a finite system replicated by periodic boundary conditions and that more computational effort is necessary compared to cutoffs (the straightforward implementation scales as N 3/2). The latter deterred most biomolecular simulators from using the Ewald method until recently, when the particle mesh Ewald method (PME) [44] was introduced. In PME, the reciprocal space part of the Ewald sum is evaluated using fast Fourier transforms, resulting in an overall N log N scaling. The potential deleterious effects of the strict periodicity that results when lattice sums are used have not been assessed in membrane simulations. However, it has been explicitly demonstrated that the use of spherical truncation introduces serious artifacts into a variety of structural, thermodynamic, and transport properties of interfacial systems and lipid bilayers [45], even when large cutoff radii are used.