膜蛋白的分子动力学研究经典文献

Setting up and running molecular dynamics simulationsof membrane proteinsChristian Kandt, Walter L. Ash, D. Peter Tieleman¤Department of Biological Sciences, University of Calgary, 2500 University Drive NW, Calgary AB, Canada T2N 1N4Accepted 17 August 2006AbstractMolecular dynamics simulations have become a popular and powerful technique to study lipids and membrane proteins. We present some general questions and issues that should be considered prior to embarking on molecular dynamics simulation studies of membrane proteins and review common simulation methods. We suggest a practical approach to setting up and running simulations of membrane proteins, and introduce two new (related) methods to embed a protein in a lipid bilayer. Both methods rely on placing lipids and the pro-tein(s) on a widely spaced grid and then ‘shrinking’ the grid until the bilayer with the protein has the desired density, with lipids neatly packed around the protein. When starting from a grid based on a single lipid structure, or several potentially di V erent lipid structures (method 1), the bilayer will start well-packed but requires more equilibration. When starting from a pre-equilibrated bilayer, either pure or mixed, most of the structure of the bilayer stays intact, reducing equilibration time (method 2). The main advantages of these methods are that they minimize equilibration time and can be almost completely automated, nearly eliminating one time consuming step in MD simulations of membrane proteins.© 2006 Elsevier Inc. All rights reserved.Keywords:Molecular dynamics simulation; Membrane proteins; Molecular modeling; Electrostatics; Algorithms1. Introduction1.1. Membrane proteinsA fundamental precondition for life is the compartmen-talization of cells and organelles from their environment that is provided by biological membranes. Beyond this basic role as a barrier, biomembranes facilitate a number of other functions that are mainly determined by the type of proteins associated with or embedded in the bilayer. Mem-brane proteins are therefore of high biological relevance, as they are key players in crucial processes such as energy con-version, transport, signal recognition, and transduction.Compared to soluble proteins it is particularly di Y cult to obtain high-resolution structural experimental data about membrane proteins: to date the 3D structures of 97 di V erent membrane proteins have been solved, whereas about 13,000 structures are available for soluble proteins [1]. Gaining a deeper insight into the architecture of mem-brane proteins is thus clearly a major goal in modern struc-tural biology. Apart from further development and improvement of structure resolving experimental tech-niques, theoretical and computational methods have become increasingly important. These range from bioinfor-matics and homology modeling procedures [2–4] to quan-tum mechanical calculations and molecular mechanical simulations [5–12].1.2. Molecular dynamics simulationsMolecular dynamics (MD) simulations numerically investigate the motions of a system of discrete particles under the in X uence of internal and external forces. The spectra of possible applications based on this approach is*Corresponding author. Fax: +1 403 289 9311.E-mail address:tieleman@ucalgary.ca (D. Peter Tieleman).476 C. Kandt et al. / Methods 41 (2007) 475–488broad, ranging from atoms in a molecule to stars in a gal-axy [13]. However, the underlying principles are the same:interactions of the respective particles are empirically described by a potential energy function from which theforces that act on each particle are derived. With knowledgeof these forces it is possible to calculate the dynamic behav-ior of the system using classical equations of motion, in their simplest form Newton’s law, for all atoms in the sys-tem. For biomolecular systems, a discrete time step of up toa few femtoseconds is used, with typical simulationsconsisting of millions of steps.For an atomic system, the potential energy function con-sists of a set of equations that empirically describe bondedand non-bonded interactions between atoms. This energyfunction together with the set of its empirical parameters is referred to as the “force W eld.” Molecular dynamics force W elds usually consist of two major components. The W rst part describes interactions between atoms connected viacovalent bonds, which typically includes bonds, bond angles, and dihedrals. The second part treats non-bonded interactions, typically as electrostatic interactions between the (partial) charges on each atom and a Lennard-Jones potential to model dispersive van der Waals interactions. Each molecule in the simulation is described by its ‘topol-ogy’, the combination of the set of all atoms with their non-bonded and bonded interaction parameters and the connectivity of those atoms in the molecule.Since their introduction in the late 1950s [14,15] andtheir W rst application to a protein [16], MD simulationshave become a common tool to investigate structure–activ-ity relationships in biological macromolecules. Providing the element of dynamics at an atomic level of detail, these simulations facilitate the interpretation of experimental data and give access to information not directly accessible by experiments. For biomolecular MD simulations, systems of more than 100,000 atoms simulated on a time scale of nanoseconds have become standard. MD simulations also form a critical component of atomic-resolution structure determination methods such as X-ray crystallography and NMR. F or a more detailed introduction to MD simula-tions see for example [13,17–19].1.3. MD simulation of membrane proteinsSince the W rst simulation of a membrane-embedded pep-tide in 1994 [20] and an integral membrane protein in 1995 [21] molecular dynamics simulations of membrane proteins have come a long way. To date, successfully investigated systems include interface-associated and trans-membrane peptides, fusion proteins, channel and pore proteins, trans-porters, ion pumps, ATP-synthases and G-protein coupled receptors. Recent overviews of the W eld can be found in [22,23].In this paper we review the most commonly used meth-ods to set up molecular dynamics simulations of membrane proteins and suggest a practical approach to setting up and running such simulations. We introduce two new and e Y cient techniques to generate a starting structure and also discuss some general questions and issues that should be considered prior to any modeling or simulation. In our group we mainly use GROMACS [24,25], a popular, freely available, and relatively user-friendly set of programs for MD simulations. At the time of writing this paper, the cur-rent version is 3.3.1. Although some of the details will be di V erent for other software packages, and for other ver-sions of GROMACS, the steps and algorithms we describe will be similar and can be readily implemented.2. Description of methods2.1. Simulation at all?When deciding to approach a speci W c problem in mem-brane protein biology using simulations, it is worthwhile to re X ect on some primary questions prior to any modeling or simulation. Is the problem to be investigated accessible by MD at all? Are the simulations practical or likely to deliver interesting and informative results? F or large membrane protein systems with a total size of 100,000 atoms or more, simulation times on the scale of tens of nanoseconds are currently accessible using advanced supercomputers. Is this time scale enough for the particular system of interest? On what time scales do the phenomena of interest take place? Are current simulations accurate enough? Small di V erences in binding free energies of the order of 1 kT are not likely to be accurately reproduced in simulations, even if there are no issues with insu Y cient sampling. Is there further experi-mental data available that can be used for validating simu-lation results? At times, the answer to these questions should suggest that simulation is not an appropriate method to address the problem at hand.If simulation seems an appropriate approach, additional questions arise. Are appropriate starting structures avail-able? Do computational models of the lipids of interest exist, or do they have to be developed (useful, but a sub-stantial e V ort)? Are unbiased equilibrium MD simulations an appropriate choice or would non-equilibrium techniques be more suitable? For example, large scale conformational changes can be forced to take place on an accessible time scale using controlled MD techniques [26,27] or by per-forming pulling experiments [28], although both methods also introduce artifacts that need to be evaluated. F ree energy-related properties can rarely be determined from equilibrium simulations of complex systems, but specialized methods, such as umbrella sampling, can sometimes be used [29,30]. Is there a parameter set that adequately describes the lipids and protein of interest, and perhaps small-molecule ligands? Many membrane proteins are sen-sitive to the lipid environment, and thus may be critically a V ected by the membrane model used [31]. Are there parameters missing that have to be determined W rst? Are topologies available for all of the system’s components? Which software should be used? This depends on the prob-lem at hand, the availability of local expertise, computa-C. Kandt et al. / Methods 41 (2007) 475–488477tional facilities, and other factors. Is enough computer time available for the planned simulations? Below we address some of the technical questions, but there is no general answer for the suitability of simulations for a particular problem.2.2. Force W elds2.2.1. OverviewForce W elds could be considered the primary assumption in MD simulations—they describe the interaction between atoms. There are many di V erent force W elds [32], but cur-rently there are only four widely used force W elds for simu-lating biological macromolecules: Amber [33,34], CHARMM [35,36], GROMOS [37] and OPLS [38]. Each of these has undergone continuous development as parame-terization methodology and experimental techniques advance, and to correctly specify a force W eld an exact revi-sion number is required. All four force W elds reproduce many protein characteristics satisfyingly well [32,39]. They have some di V erent strengths and weaknesses, but most importantly they share a number of the same limitations due to the basic simplifying assumptions necessary in large scale MD simulations. The treatment of electrostatics in particular is somewhat simplistic because electronic polar-izability is only accounted for in an average way. F orce W eld and methods development aimed at resolving these de W ciencies is ongoing [32,39–41].Phospholipids are a subset of biomolecules that have received considerably less attention than proteins or nucleic acids. There are only two phospholipid force W elds in com-mon use today, although additional sets are in development (e.g. [42]). The W rst commonly used lipid force W eld is part of CHARMM [43,44]; the second is based on an older ver-sion of OPLS and AMBER with some additional parame-ters from Berger et al. [45]. Both force W elds are able to reproduce the available information from stable bilayers fairly well, although both sets of lipid parameters will ben-e W t from further re W nements. Recent examples are the treat-ment of the rotational isomerization of acyl tails, where intermediate-range intramolecular interactions play an important role [46], and a direct comparison of simulations of lipids with di V raction experiments [47].Unlike CHARMM lipids and most of the recent protein force W elds, Berger lipids are described by a united-atom force W eld; each aliphatic carbon with associated hydrogens is described by a single particle with the approximate physi-cal characteristics of a methyl, methylene, or methine group. This is a level of detail somewhere between fully atomistic descriptions and some of the coarse-grained approaches that have been developed over the years [48,49]. When hydrogens are treated explicitly, the number of atoms per lipid approximately triples, and the number of pairwise interactions in the membrane becomes much higher. Because the computationally inexpensive united-atom lipids reproduce experimental behavior reasonably well, it may be desirable to use a combination of united-atom lipids and an accurate all-atom protein force W eld [50].Ultimately, an all-atom force W eld should be more accurate,but despite some discussion in the literature we do not believe that we have reached the point yet where the united-atom approximation is the limiting factor in accuracy formost purposes.This point is related to a broader problem. Signi W cant progress has been made in developing modern protein force W elds, based on a comparison with high-level quantum mechanics and high-resolution experiments on soluble pro-teins, but it is challenging to obtain accurate experimental data to validate and improve lipid models. The situation is even more complex for lipid–protein interactions compared to pure lipids, because experiments on membrane proteins and peptides typically have a signi W cantly lower resolution than in solution. In addition, from a simulation point of view, critical tests are di Y cult due to the long simulations required. F or now, the most straightforward approach to membrane protein simulations involves directly combining mathematically compatible protein and lipid force W elds, and this approach has yielded useful insights into membrane protein behavior [22,23].A key di Y culty with classical simulations of molecules that move between two very di V erent environments is that the parameters used to describe them have to be accurate in both environments, although typical biomolecular force W elds have been developed for aqueous solution. A direct approach to testing the free energies of transfer between water and hydrophobic environments is now possible com-putationally [50–54], which allows a direct comparison to experimentally measured partition coe Y cients [55] and the interfacial partitioning of model peptides [56]. Such studies have found substantial errors in free energies of transfer for amino acid side chain analogues between water and cyclo-hexane. Signi W cant improvements in current force W elds can be expected based on re-parameterization e V orts, both in the thermodynamics of membrane proteins at the water/ lipid interface and in processes such as small molecules binding to proteins, which typically also involves a change of environment. There will be a limit to the maximum accu-racy that can be reached by simple re-parameterization due to ignoring electronic polarizability e V ects (e.g. [40]), but in our opinion we have not reached the point yet where fur-ther improvement is impossible without including such e V ects. An intriguing question is whether a general polariz-able force W eld would solve many problems at once, possi-bly with ultimately less e V ort. Undoubtedly we will see many endeavors in the near future aimed at improving protein and lipid force W elds.2.3. Building topologiesMembrane proteins are often associated with speci W c cofactors, ligands or other biological molecules for which the simulation package holds no topologies. In many cases, such descriptions might already exist in another force W eld and might be converted; in other cases, they must be built478 C. Kandt et al. / Methods 41 (2007) 475–488anew. However, care must be taken to ensure their accuracy and compatibility with the protein and lipid models. A use-ful and commonly used approach when generating a new topology is to dissect the molecule into more manageable parts. For example, complex substrates and cofactors often contain amino acid-like fragments for which force W eld parameters are already available, which thus require only slight adjustments. High-resolution X-ray data can be used for reference when parameterizing bond length and angle de W nitions. F orce constants for bonds and angles can in principle be obtained from accurate QM calculations, but in practice bonds are often taken as W xed (so that only the bond length matters) and angle force constants are chosen similar to existing angle parameters in the force W eld.One important factor to consider when evaluating a model for a cofactor is its partial charge distribution, or more generally, the way the force W eld models non-bonded interactions. Point charges and Leonard-Jones interactions that give the proper behavior for a molecule in water may not be adequate to describe behavior in the interior of a bilayer or a hydrophobic pocket in a protein. The charge distribution (usually derived from ab initio QM methods) should be compatible with the protein force W eld; a charge partitioning technique that overestimates charges on a cofactor relative to a protein may arti W cially and adversely a V ect the balance of forces controlling protein-cofactor interactions. The interactions between ions and a potassium channel are an interesting example [57]. Thus it is advisable to use a methodology similar to that used in parameterizing the protein (and/or lipid) force W eld.Topologies of complex molecules present a signi W cant practical problem. MD software typically generates topolo-gies for biopolymers like proteins and polynucleotides automatically from the 3D coordinates, but this relies on an existing description (in the software) of the exact topology of amino acids and individual nucleotides. It is a challeng-ing problem to generate accurate topologies for arbitrary molecules. Although there are a number of web servers [58–60] and computer programs (e.g., AMBER, Insight II) that will attempt to generate topologies based on chemical structures or 3D coordinates, in practice these require care-ful manual checking and sometimes adjustments. At the moment topologies of complex molecules like vitamin B12 require signi W cant manual e V ort [61].2.4. Creating a starting structureIn addition to addressing general questions and chosing simulations parameters, setting up a membrane protein simulation system usually requires four major working steps, which will be discussed in the following sections. These steps are: 1, orienting the protein; 2, preparation of the bilayer; 3, solvating the system; 4, system equilibration. As pointed out earlier, in describing these processes we make explicit reference to GROMACS utilities [24,25] dis-tributed at the time of writing (version 3.3.1). However the methods described herein do not rely on speci W c software and can be readily implemented with other molecular dynamics and modeling programs.2.4.1. Orienting the proteinWhen starting from a pre-existing bilayer the common approach begins with centering the protein in the simula-tion box. Both the protein and the bilayer box should have the same dimensions so the two coordinate sets can be com-bined easily. The protein is then oriented in such a way that its hydrophobic belt is aligned with the non-polar lipid tails. This can be done using either common molecular graphic programs such as MolMol [62], VMD [63] or RasMol [64] or via the GROMACS tool editconf, applying explicit translation and rotational vectors to the protein only. The term hydrophobic belt refers to a common feature of mem-brane proteins regarding the distribution of charged resi-dues on the molecule’s surface: the area exposed to the hydrophobic component of a bilayer usually contains no charged residues [65,66] (see Fig.1).For peptides the procedure is the same or even simpler when dealing with short (»20 residues) trans-membrane peptides, where their entire length does not exceed the bilayer thickness. When constructing the membrane around the protein instead of starting from a pre-equilibrated bilayer, one should ensure that lipids are placed in such a way that the hydrophobic belt criterion is met. F urther re W nements are possible, as the hydrophobic belt criterion is not exact and could be replaced by a more quantitative assessment based on e.g., optimizing the electrostatic free energy of solvation using continuum electrostatics [67].2.4.2. Preparation of the bilayerMany important motions exhibited by lipids and lipid bilayers are relatively slow compared to the time scale of MD simulations [68–70], including such processes as lipid di V usion and rotational reorientation, phase transitions, bilayer self-assembly, and lipid X ip-X op. Therefore, one major di Y culty in designing a membrane protein simulation Fig.1. (a) Membrane proteins typically exhibit a characteristic distribu-tion of charged residues (black) on their molecular surface (white) whereas the area exposed to the hydrophobic component of the bilayer consists of neutral residues only (a). This hydrophobic belt (b) can be used to align the protein with the polar lipid head groups (grey). (F igure created using MolMol and XaraX1).C. Kandt et al. / Methods 41 (2007) 475–488479experiment is the proper assembly of the starting lipid bilayer and protein system. There are, in essence, two possi-ble approaches to this problem: 1. A protein is inserted into an existing, equilibrated bilayer. 2. A bilayer is built around the protein to embed it in the bilayer.We W rst review the main methods used in the literature. Then we introduce two new and e Y cient approaches that we recommend for new simulations: one starting from the coordinates of a single lipid molecule and one starting from a pre-equilibrated bilayer. Either way, if not stated other-wise all preparation techniques initially take place in the absence of water.In the end it does not matter which method one uses to make a starting structure as long as the system has been equilibrated long enough and known experimental data are adequately reproduced. However, the practical di V erences between approaches can be substantial. Important criteria to favor one method over another are the time required to equilibrate the system, the amount of manual e V ort required, and the potential for serious artifacts due to deformations of the bilayer caused by the wrong number of lipids in one of the lea X ets relative to the other lea X et.2.4.2.1. Review of preparation techniques2.4.2.1.1. Delete lipids within a cut-o V. This is the sim-plest and most commonly used method of preparing a bilayer for membrane protein simulation (Fig.2a). Recent examples can be found in [71–80]. The coordinate sets of bilayer and protein have already been combined and over-lapping lipids are removed on the basis of a simple distance cut-o V between the protein and either single lipid atoms (usually the phosphorus in the head groups), or the entire lipid. The applied cut-o V length varies between 5–6Å (based on protein–phosphorus distances) and 0.8–1.6Å (based on any contact between protein atoms and whole-lipids).The main problem with this method is that, due to the highly disordered lipid conformations, a rather rugged hole is created which can require a long simulation time before the local density has returned to its equilibrium value again. This can also lead to a shrinking of the whole membrane patch below an area that is not suitable for the system size anymore. The graphical representation option “draw peri-odic images” in VMD o V ers a particularly useful way to check if the chosen system size is big enough to ensure the protein does not “see” any of its periodic images in the course of the simulation. It is generally useful to keep in mind that lipid di V usion takes place on a time scale of 10–100ns and for typical simulation studies equilibration times of 10–20ns are necessary [68]. Non-equilibrium situations like the presence of a hole might initially decay rapidly, but the resulting structures may still not represent equilibrium for a signi W cant amount of time.Another problem frequently encountered, especially when using the protein – phosphorus approach, is the potential for clashes between protein residues and those lip-ids which are still present after the removal step. One possi-ble way to deal with this is to exchange the a V ected protein side chains with alternative and more suitable rotamer con-formations. This can be done using packages such as Swis-sPDB [81,82], but we prefer to avoid this situation altogether.A variant of the cut-o V method uses a rectangular grid placed on bilayer and protein. All lipids with atoms fall-ing in grid cells already occupied by protein are removed [83].2.4.2.1.2. Using repulsive forces to create a hole – cylinder & mdrun_hole. The protein shape is approximated by a cyl-inder of appropriate radius and all lipids with their phos-phorus atoms inside this cylinder are removed (Fig.2b). Subsequently all lipid molecules which still have (non-phos-phorus) atoms inside are driven out by applying a weak repulsive potential, after which the protein can be inserted in the generated cavity [84]. The major drawback of this method is clearly the large simpli W cation of protein shape. While it may work W ne for single helices or highly symmet-rical helix bundles [85,86], it is not suitable to deal with membrane proteins of more irregular shaped cross sections. The program code has been implemented in the GRO-MACS package and recent applications of this method can be found in [87–89]. We consider the newer mdrun_hole (see below) procedure superior.Introduced in 2002 [90], mdrun_hole is a further devel-opment of the cylinder approach. One still starts with a cyl-inder as a W rst approximation of the protein shape and deletes all lipids with phosphorus inside that cylinder, but then the Connolly surface of the protein itself is used as ref-erence boundary to drive out lipids and water molecules (see Fig.2c), creating a tailored cavity to contain the pro-tein [91–98]. The program code is part of the GROMACS package and additionally requires GRASP [99] or another program to generate the solvent-accessible protein surface in a format that can be read by GROMACS. The process leads to systems that have excellent lipid–protein packing and require minimal time to equilibrate compared to coarser methods. Cylinder and mdrun_hole can be applied to both water-free and fully hydrated systems.In practice, a drawback of this procedure is that manual intervention is often required: mdrun_hole has quite a few input parameters that need to be speci W ed but there is no standard setting that works for every system. Instead parameters like force constants for repulsive potential, sur-face de W nition (atoms to consider, probe radius), o V sets (de W ning the e V ective position of the surface boundary), and water repulsion parameters may need to be W ne-tuned for each system individually. Lipids can also become trapped on their way out, in small cavities created by clus-ters of side chain atoms, and must be removed manually. If parameters for pressure and force constants are not chosen properly (constant pressure and appropriately strong force constants) and/or the number of lipids deleted initially is too small, it is possible that some of the lipids may X ip over, with their tails pointing towards the water phase. Con-straints on lipid z coordinates can be used to prevent lipid X ipping.480 C. Kandt et al. / Methods 41 (2007) 475–488Fig.2. As discussed in the text the preparation of a bilayer can be done by deleting lipids within the range of a simple protein – lipid distance cut-o V (a), by applying a repulsive potential to drive out lipid and water from the area to be occupied by the protein as in the cylinder approach (b), or with mdrun_hole (c). Alternatively the membrane can be constructed around the protein lipid by lipid, taking structures from a pool of hydrated lipid conformations (d). We introduce two approaches that involve scaling of the simulation system: one method starts by placing copies of lipids on a 2D grid and then scales down the box size stepwise with energy minimization and short molecular dynamics steps in between (e). The other uses a pre-equilibrated bilayer, which is expanded within the xy-plane. Lipids within a given cut-o V distance are removed and the system is brought back to natural dimension by a series of scal-ing steps of the lipid xy positions and energy minimizations (f). (Figure created using RasMol, VMD and XaraX1).。