Bioinformatics for Membrane Lipid Simulations: Models, Computational Methods, and Web Server Tools

Abstract

Biological membranes are complex environments consisting of different types of lipids and membrane proteins. The structure of a lipid bilayer is typically difficult to study because the membrane liquid crystalline state is made up of multiple disordered lipid molecules. This complicates the description of the lipid membrane properties by the conformation of any single lipid molecule. Molecular dynamics (MD) simulations have been used extensively to investigate properties of membrane lipids, lipid vesicles, and membrane protein systems. All-atom membrane models can elucidate detailed contacts between membrane proteins and its surrounding lipids, while united-atom and coarse-grained description have allowed larger models and longer timescales up to microsecond mark to be probed. Additionally, membrane models with mixed phospholipids and lipopolysaccharide content have made it possible to model improved views of biological membranes. Here, we present an overview of commonly used lipid force fields by the biosimulation community, useful tools for membrane MD simulations, and recent advances in membrane simulations.

1. Introduction

The biological membrane is a selective barrier delineating the boundaries of cells and organizing cellular organelles into compartments. One of the major functions of the membrane is to regulate movements of water, ions, and certain constituents in and out of the cell or organelle. Singer and Nicholson described the membrane as a “fluid mosaic” model in which lipids and proteins freely diffuse in the membrane plane [1]. While the biological membrane is actually made up of 1.5- to 4-fold proteins by weight, it is usually described as phospholipids arranged in a bilayer.

The lipid composition in membrane is highly specific with variation identified from membranes of different organisms, different cells of the same organism, and even different membranes of the same cell [2]. Generally, a lipid consists of a polar head group and hydrophobic tail region where a cell may have over a thousand different lipid species. Formation of bilayers and other lipid structures such as micelles occurs as the hydrophobic tails orientates away from the aqueous environment while the polar head groups interact favorably with water. Changes in the bilayer structures can affect the function of the membrane, membrane-embedded proteins and complexes. The surrounding environment can greatly affect the structure and functions of proteins, much like the effects of water on water-soluble proteins. The membrane produces a complex environment for embedded proteins, and membrane thickness, fluidity, charge, curvature, and phase have been demonstrated to play a role in protein structure and function determination [3].

The structure and properties of membranes are a complex matter and cannot be easily described by a single-lipid molecule alone. Molecular dynamics (MD) simulations offer a viable alternative to study the properties of membrane and lipid-forming structures such as vesicles. This review serves as an introduction to currently available membrane lipid force fields and recent advances in membrane simulations.

2. Force fields for lipid simulation

In general, all-atom (AT), united-atom (UA), and coarse-grained (CG) are the three-membrane lipid force fields. Figure 1 illustrates the AT, UA and CG force field of a lipid by spherical representation.

Figure 1.

Representation of 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) with (a) atomistic (all-atom; AT), (b) united-atom (UA), and (c) coarse-grain (CG) force fields as van der Waals spheres.

2.1. All-atom (AT) MD simulation: details and details

AT MD simulation represents every atom in the system as a single interaction site. Figure 2 shows the example of AT simulation system of Salmonella enterica ser. Typhi TolC protein in POPE. To date, Chemistry at HARvard Macromolecular Mechanics (CHARMM) and Assisted Model Building with Energy Refinement (AMBER) are the only fully AT force field parameterization available for lipids. Prior to the development of CHARMM36, CHARMM27r was widely used for membrane simulations [4, 5]. Simulations using CHARMM27r require a large positive surface tension (30–40 dyn/cm) to achieve the experimentally determined surface area per lipid (APL). However, theoretical considerations of self-assembled systems and macroscopic black lipid bilayers [6] indicate that the surface tension of bilayers are about zero to several dyn/cm even when undulations are taken into account [7]. Therefore, simulations of lipid bilayers using CHARMM27 and CHARMM27r shrinks to a near gel phase state without the use of surface tension. Additionally, CHARMM27 and CHARMM27r also failed to reproduce the experimental deuterium order parameters, S_CD in the glycerol and upper chain regions [4]. A wide range of glycerophospolipids exhibit splitting in the carbon 2 of the aliphatic chain and carbon 1 for glycerols, but this observation cannot be replicated when simulations were performed with CHARMM27 or CHARMM27r [8]. This may affect conclusions of interactions of lipids with surface-active agents drawn from MD simulations. Besides that, the area compressibility modulus, K_A, was underestimated, the head group region was underhydrated, the electron density in the bilayer midplane was underestimated while the frequency dependence of the ¹³C NMR T₁ of the acyl chains near the head group was overestimated.

Figure 2.

All-atom (AT) simulation of Salmonella enterica ser. Typhi TolC protein in POPE.

These major weaknesses of the CHARMM27r force field eventually motivated the development of CHARMM36. CHARMM36 corrected most of the existing problems with lipid force fields, most importantly allowing the simulation of lipid bilayer without the use of surface tension and improved the reproducibility of the experimental deuterium order parameters in the glycerol and upper chain regions of phospholipids [8]. A comparative study of lipid force fields by Piggot and colleagues found that CHARMM36 was the only force field that accurately reproduced the experimental order parameters of carbon 2 in both acyl chains of 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) and 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (POPC) [9]. The other properties of the membrane were reasonably well reproduced. However, there are several limitations that have to be considered while pursuing CHARMM36. One is the approximate treatment of long-range Lennard-Jones (LJ) forces for bilayers. Simulations augmented with long-range LJ forces using 3D-isotropic periodic sum/discrete fast Fourier transform (3D-IPS/DFFT) may improve results for lipid monolayers but this increases surface tension, therefore making it less suited for bilayer simulations. The recommendation for bilayer simulations is to use particle-mesh Ewald (PME) with r_c = 10 or 12 Å and no long-range correction for the LJ term, but these setting will cause underestimation of the surface tension of lipid monolayers.

AMBER force fields were generally less used for membrane protein simulation due to the lack of a specific parameter set for lipids. However, the general AMBER force field (GAFF) which was originally parameterized for the simulation of arbitrary organic molecules with pre-existing AMBER force fields have been shown to reproduce lipid bilayer parameters satisfactorily [10]. GAFF has been tested on a range of lipids, including 1,2-dimyristoyl-d54-sn-glycero-3-phosphocholine (DMPC), 2-dioleoyl-sn-glycero-3-phosphocholine (DOPC), 1,2-dilauroyl-sn-glycero-3-phosphocholine (DLPC), DPPC, POPC, and 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine (POPE) [11, 12]. The use of surface tension is necessary to achieve the correct APL for POPC [10], DOPC [13], and DMPC [12]. In order to better fit GAFF for lipid molecules, Dickson and colleagues set to re-parameterize the LJ terms for acyl chain carbons and hydrogens as the initial GAFF LJ terms were developed for proteins, nucleic acids, and small organic molecules [11].

In addition, torsion parameters were also re-optimized using high-level quantum chemical data with “Paramfit.” This strategy has allowed a tensionless simulation of the lipids in the isothermal-isobaric (NPT) ensemble, achieving high level of agreement with experiment data for volume per lipid (V_L) and thickness value within 5% of experiment [11]. However, the APL for POPE lipids is lower than expected and thickness value was overestimated by 10% [11]. The force field has also been shown to be able to reproduce the data on large lipid bilayers containing 288 lipids, with little changes in the APL, V_L, and bilayer thickness [11]. The comparison with CHARMM27 and Berger force field also displayed GAFF's ability to reproduce the order asymmetry found at the beginning of lipid acyl chains [13]. Prior to the development of CHARMM36, GAFF was the only force field capable of capturing the carbon-2 deuterium order parameter splitting [13], although it has also produced lower surface APL and higher deuterium order parameters compared to experimental data [10, 12, 13]. Even so, compatibility of GAFF with previous AMBER force fields makes it a suitable choice for the simultaneous simulation of membrane, protein, and organic molecule simulation [13].

A recently developed Lipid14 force field with an updated Lipid11 head group and tail group charges and parameters enabled proper tensionless simulation of lipid bilayers in AMBER [14]. Lipid14 LJ and torsion parameters were modified to reproduce the experimental density, ρ and heat of vaporization, ∆H_vap of alkanes of different chain lengths by fitting the CH₂-CH₂-CH₂-CH₂ torsion to ab initio data using “Paramfit” and altering the LJ and torsion simultaneously. Lipid14 force field has a modular nature that allows new lipid species to be added into the force field by constructing them from head group and tail group “building blocks.” Testing of Lipid14 on DLPC, DMPC, DPPC, DOPC, POPC, and POPE showed that APL for all simulations was within 3% of the experimental value, with the exception of POPE. The APL was closer to the older experimental value of 56.6 Ų [15] than the more current APL of 59–60 Ų [16]. The V_L was found to be within 5% of experimental value which was acceptable, although it may be considered a slight underestimation [14].

Furthermore, the isothermal area compressibility modulus, K_A falls close to the experiment, again with the exception of POPE, which has a higher K_A with a large standard deviation value. The authors suggested that implementation of other barostats into AMBER may improve K_A values as the Berendsen method for pressure control is not ideal for simulations in which volume fluctuations is an important parameter that is capable of influencing the outcome of the simulation [14]. The Luzzati thickness, D_B, which is calculated using the z-dimension of the simulation box and the integral probability distribution of the water density along the z-axis was slightly underestimated for saturated lipids, which implies more water penetration into the hydrophobic region of the bilayer of these lipids. Lipid14 was also able to reproduce the experimental order parameter trend, including the splitting of the sn-2 chain from the sn-1 chain and the drop at the carbon-9 and carbon-10 positions of POPC and POPE lipids, which is the cis double bond. Results from GPU repeats and CPU runs were also consistent.

Recent analysis applying full AT force field have reached timescales of up to tens of nanoseconds, with ambitious simulations pushing the microsecond mark, e.g. to characterize the interaction of the multiple sclerosis synthetic biomarker CSF114(Glc) with the membrane bilayer [17], probing of the huntingtin Htt17 membrane anchor on a POPC bilayer [18], mixing of lipids [19], and membrane-binding mechanism of the yeast Osh4 peripheral membrane protein [20]. Even though microsecond simulations have been reported, the accessibility of this approach is limited as many biologically relevant phenomena may require sampling time beyond the microsecond timescale and investigators may not have access to intensive computational resources. Given that lipid reorganization is quite slow, a properly equilibrated membrane may require 20–40 ns of simulation time. Convergence of membrane protein simulations still remains a question, as 100 ns of simulation may still be insufficient to fully describe e.g. rhodopsin loop fluctuations in a membrane [21]. Indeed, since all atomic-level interactions are retained and time-steps for integration of Newtonian motions are in the femtosecond range, AT simulations could be a time consuming and computationally expensive practice. The issue of convergence has been rectified in part by running several long simulations and using non-equilibrium sampling method or steered MD to calculate the observables of interest. A series of extended simulations (~100 ns) on the voltage sensor domain of potassium channels revealed the importance of lipid phosphates in accommodating the significantly charged S4 helix [22]. Using umbrella sampling, potential of mean force has been calculated for permeation and effect of a potassium ion on the second ion passing through gramicidin A channel [23].

Regardless, AT simulations still provided the highest level of detail and reliability when it comes to quantitative prediction of properties such as motional timescales or interaction strengths. Such details recently provided insightful interactions between protein and membrane, such as the rearrangement of amino acid side chains and local bilayer deformation due to hydrophobic mismatch and how the hydrophobic membrane layer accommodates charged arginine side chain of outer membrane phospholipase A [24, 25].

2.2. United-atom (UA) lipid models: the best of both worlds?

The UA representation of lipids simplifies the carbon tails of the lipid by associating the aliphatic carbon and its hydrogen atoms into a single particle. Because the non-polar hydrogen atoms are treated implicitly, the number of interaction sites per lipid can be reduced by two third. The computational costs for simulations of such membrane systems become relatively cheap as the 60% of the pairwise interactions in the membrane is reduced. The model lipid DPPC can be represented by 50 particles in UA force field, but needed 130 interaction sites in an AT force field. Since limited physical information may be collected from explicit acyl chain hydrogens, it became desirable to utilize UA force fields for membrane with an AT protein force field for membrane protein simulation [26, 27].

The UA lipid models parameterized by Berger et al. (1997) were one of the most popular lipid force field for lipids and were originally developed by Essex and colleagues [28] from the Optimized Potentials for Liquid Simulations (OPLSs) UA force field. Bonded parameters of the Berger lipids were obtained from the GROMOS87 force field (note: GROMOS is the GROningen Molecular Simulation package), the acyl chains used Ryckaert-Bellemans dihedral parameters whereas the van der Waals terms were from OPLS and atomic partial charges were from Chiu and colleagues' calculations [29]. Berger and colleagues further optimized the LJ parameters of the lipid tails based on thermodynamic data of pentadecane [30]. Berger lipid parameters were recommended if one desires maximal sampling due to its fast diffusion and good simulation efficiency [9]. For membrane protein simulations, Berger lipids are commonly used with OPLS and GROMOS [27]. It has also been demonstrated to be compatible with AMBER99SB and could give marginally better free energy calculations than the widely used OPLS/Berger combination [31]. While simulations of membrane proteins using such hybrid parameters have been validated by various groups [32–34], the combination of different force fields require care. Protein lipid interactions in Berger-OPLS and Berger-GROMOS have been found to be overestimated and result in drastic changes of lipid properties upon protein insertion [13, 27].

A CHARMM UA representation of the acyl chains is also available, and compatible for simulations with AT CHARMM protein force fields [26]. The CHARMM UA model was derived from C27 phospholipids, where explicit hydrogen atoms of the acyl chains were replaced with a UA representation. The force field, called C27-UA, retains the accuracy of the AT counterpart and provides a practical alternative when used in simulations of proteins and other compounds described by C27. C27-UA was parameterized by fitting to experimental data and AT simulations of liquid phase model systems (pentadecane for saturated hydrocarbon chains, cis-5-decene for monounsaturated chains and methyl hexanoate for the ester region). Simulation of POPC bilayer with C27-UA was comparable to C27 and reproduced experimental NMR and X-ray diffraction data, including electron density profile and carbon-deuterium order parameter. The free energy profiles of transfer of ethane, methanol, and water across a water-dodecane interface were identical in C27-UA and C27 simulation, suggesting that the force field is capable of simulating mixed system containing UA lipids and AT proteins and organic molecules described by the standard CHARMM force field. However, it also retained C27's feature of requiring a positive surface tension to be able to reproduce the experimental APL.

Several GROMOS parameters are used in membrane simulations, such as 43A1-S3 [35], 53A6 [36], and Kukol's modification of the 53A6 parameters [37]. 43A1-S3 is an extension and modification of the 43A1 force field designed to improve the properties of lipid membranes in simulations. 43A1-S3 employs charges from Chiu et al. [29] while van der Waals and dihedral parameters were modified from 43A1 to improve hydrocarbon and choline head group dynamics [38]. GROMOS 43A1-S3 accurately reproduces many properties of DPPC bilayers, including APL, lipid diffusion, and the order parameter [39]. The 53A6 parametrization also uses Chiu et al.'s charges [29] and LJ parameters of the choline methyls and phosphate ester oxygen atoms, thus providing good agreement for the APL, density profile, and order parameter for saturated and unsaturated acyl chain PC lipids [40, 41]. Even with improvements introduced, the isothermal area compressibility modulus was still overestimated. Meanwhile, Kukol reparametrize 53A6 following reports that the force field failed to reproduce DPPC satisfactorily [37]. He used Chiu's charges and increased the van der Waals radii of the carbonyl carbons of DPPC, DMPC, POPC, and POPG. He also employed non-standard GROMOS dihedral parameters for the double bond in the unsaturated lipids. This results in fairly good agreement with experimentally determined properties with the exception of order parameters for the sn-2 oleoyl chains. A comparative force field study by Piggot et al. recommended against using 43A1-S3 for POPC membranes and Kukol's 53A6 POPC parameters [9]. This is because 43A1-S3 could not reliably reproduce the drop in order parameter at the carbon 10 of oleoyl chains whereas Kukol's POPC parameters showed several disagreements with experimental value in terms of membrane thickness and order parameter [9].

With the advances in computational power, large-scale simulation projects using UA representation of the membrane region have been able to reach the microsecond mark. Access to such timescale allowed the observation of structural rearrangement and changes in hydrogen bonding pattern in an integral Kv1.2 channel embedded in a hyperpolarized membrane [42]. On a shorter timescale, CorA magnesium transport channel has been probed to undergo conformational changes to a putative open state in 110 ns [43].

2.3. Coarse-grained (CG): the need for speed

CG simulations are being widely used to investigate phenomenon occurring in timescales not accessible by AT simulation. In a CG simulation, 3–4 heavy atoms (non-H) are grouped together and represented by a single particle. For example, a DMPC lipid consisting of 130 atoms can be represented by 12 interaction sites [44]. The choline moiety is modeled with a single positively charged particle, the phosphate group with a negatively charged particle, the glycerol linkages with two nonpolar particles, while the lipid chains are modeled with 4 apolar particles each [45].

Early CG approaches were typically parameterized based on comparison to AT simulations by using inverted Monte Carlo schemes [46–48] or force matching [49], which aims to reproduce the structural details at a particular state for a particular system. Instead, the MARTINI (note: MARTINI is a CG force field developed by Marrink and coworkers [44]) CG force field calibrated the building blocks of biomolecules against thermodynamic data, particularly the oil/water partitioning coefficients so that there is no need to re-parameterize the model each time. A consistent, atomic level compatible CG force field will be beneficial for multi-scale applications. In MARTINI, an average of four heavy atoms were represented by a single interaction site, with the exception of ring structures which has 2 or 3 ring atoms mapped to a CG bead. The beads are then distinguished into four main types, i.e., polar, nonpolar, apolar, and charged. MARTINI was able to reproduce the properties of lipid bilayers on a semiquantitative level. These are including the APL, the distribution of groups across the membrane and the bending and area compression moduli.

MARTINI 2.0 further improved the stress profile across the lipid bilayer and its tendency to form pores. Additionally, the free energy of lipid desorption and lipid flip-flop across the bilayer agreed well with AT simulations and the condensing effect of cholesterol on the APL can be reproduced. MARTINI achieved 5- to 10-fold faster sampling of the configurational space of liquid hydrocarbons and the lipid tails inside a bilayer compared to AT force fields [44]. Among others, MARTINI has allowed applications in simulations of vesicle formation and fusion [50, 51], phase transition of lipid bilayers [52], and the structure and dynamics of membrane-protein assemblies [53, 54]. Currently, MARTINI has also been supplemented with an implicit solvent force field, named Dry MARTINI, which provided 1–2 order speed up and is expected to find application in simulations of large membranes containing millions of lipids [55].

The ELBA (an acronym for “electrostatic-based” by Orsi and Essex [56]) force field, so named because it is electrostatics based, offers another alternative for CG simulations of lipids and lipid bilayers. The ELBA model features two approaches in contrast to other CG methods. Notably, LJ interactions are treated using standard Lorentz-Berthelot mixing rules, similar to AT force fields. This is due to the explicit treatment of lipid electrostatics and water dipoles whereby a relative dielectric constant of unity (ɛ_r = 1) was used to model their interactions. In addition, a realistic dipolar water model based on a simpler soft sticky dipole potential, i.e. the Stockmayer potential is used in ELBA, which helped provide a correct diffusion coefficient of lipids in the liquid phase, a point often overestimated by other CG models. It also shows 15 and 200 times speed improvement over UA and AT models, respectively. Validation has been performed on DOPC, DOPE, and gel phase DSPC, whereby ELBA satisfactorily reproduced several fundamental experimental properties including APL, volume per lipid, curvature elastics constants, electrostatic potential distribution, electrostatic diffusion constant, and lipid diffusion coefficient. At the time of the development, ELBA was only available with the authors' in-house software, and it has since been implemented in LAMMPS [57]. However, there has yet to be a current parameterization for protein and other organic molecules, potentially limiting its use in mixed systems.

CG models for proteins are available and have been adopted into protein-bilayer and peptide-bilayer systems [58, 59]. The advantage in employing CG models is the improved speed and size. As CG force fields have a reduced degree of freedom and removed high-frequency motions such as hydrogen bond vibrations, integration time step could be pushed up to 20–40 fs and highly increases the accessible timescale by 100-fold. Therefore, it is possible to access phenomenon occurring at timescales not accessible to classical simulations, for example, membrane protein aggregation, demonstrated by the formation and domain-specific distribution of Ras proteins in plasma membranes [60]. CG simulation has also been used to self-assemble lipid bilayers around membrane proteins, providing information into protein positioning in the bilayer [61]. Comparison of membrane positioning for six proteins in the CG approach with available experimental data showed high qualitative similarity, even though there are discrepancies from different lipid composition used in the simulation and experiments [61].

Nevertheless, due to simplification of the representations, loss of details is inevitable. Two approaches have been employed to probe for detailed interactions while simulating in CG membrane models is by mixed AT-CG simulation and multiscaling. To elucidate small molecules, membrane peptides and proteins at high resolution, some approaches have used mixed CG-AT systems by simulating the ligand and protein atomistically in a CG membrane [62, 63]. Such mixed descriptions can be likened to quantum mechanical/molecular mechanical (QM/MM) simulations that have achieved considerable success in the past decade (see [64, 65] for reviews regarding usage of QM/MM with proteins). A major consideration in building a mixed AT-CG system is the parameterization and definition of the interactions at the AT-CG interface. As many CG force fields are derived from AT simulations, force-matching procedures can be used to derive effective pairwise CG force field from AT simulations. In AT-CG approach, all atoms simulation of the whole system is first carried out, the system will be subsequently divided into AT and CG parts whereby the most interesting part of the system will retain the atomistic details. The effective AT-CG force field is subsequently obtained by treating the AT and CG parts equally in the force matching procedure [63].

CHARMM-GUI PACE CG Builder [66] was developed for the modeling of large and complex biological system. It utilizes the mixed UA/CG where the PACE force field was used for protein in UA format and Martini CG force field for the water, ions and lipids. Thus, the number of atoms in a system can be reduced by the factor of 10. Analysis showed the PACE/MARTINI hybrid simulation has most of the proteins in the root means square deviation (RMSD) of less than 3 Å.

On the other hand, several methods have been developed for conversion of CG system to AT representation through fragment-based approach [67], simulated annealing [68], and force-matching [49]. These methods, termed the multiscale approach allows the use of CG simulations to explore membrane-lipid interactions, which, after a sufficient equilibration period is converted back to AT model simulation for detailed characterization [67].

3. Tools for membrane simulation setup and analysis

3.1. Automated setup of membrane simulation systems

The availability of many different force fields and parameters for a range of lipid molecules has made it easier to construct systems consisting of mixed lipids. Therefore, the lipid-converter tool is beneficial to easily adapt a system between force fields [69]. The lipid-converter tool can be used in command line by defining a PDB or Gromacs coordinate file and is also available as a web server [70]. The tool currently supports the Berger, GROMOS 43A1-S3, GROMOS 53A6, GROMOS 53A7, CHARMM36, OPLS-UA, and Lipid11. Thus, Stockholm lipids that are compatible with AMBER force field and used CHARMM nomenclature are also supported by this extension. Moreover, lipid converter may be useful to build non-conventional systems as it can generate asymmetric lipid distribution and even label leaflets in curved systems like vesicles.

Towards the end of automating the process of building heterogeneous membrane, the CHARMM-GUI Membrane Builder [71] is a useful tool to generate coordinates for membrane models and protein/membrane systems [72–74]. The membrane builder offers a selection of commonly used lipid models in addition to cholesterol which can be customized according to concentration, APL, hydration number, and thickness of the water layer [73].

Alternatively, a new web server MemGen is able to automatically set up lipid membrane simulation systems without restrictions in force fields, lipid types, or MD simulation software [75]. The user can upload one or more lipid structure files as well as amphiphilic molecules such as alcohol or detergent. A compact representation of each lipid aligned along the z-axis is generated by building GAFF topology of each lipid using ACPYPE, then applying simulated annealing with constant-force pulling on the head group and tail atoms as well as position-restraining potentials with Gromacs. The server subsequently hydrates the membrane with a number of water molecules, which can also be specified by the user. After the addition of counter ions or sodium chloride, a PDB format of the final structure is available for download. However, it must be noted that MemGen provides highly ordered, unphysical configurations which requires careful equilibration of at least 10 ns. In addition, it is also unable to produce asymmetric bilayers with different composition of lipids in the two monolayers.

iMembrane is another useful web-based tool which can predict the orientation of a membrane protein within the membrane [76]. Early approaches use a two-state membrane model or a simple hydrophobic slab to model the orientation of a membrane protein in the membrane. Scott and colleagues developed a CG MD to simulate membrane proteins in the presence of membrane lipids which self-assemble into a lipid bilayer [61]. Using the simulation results, iMembrane can predict the orientation of proteins of homologous structure or sequence. BLAST is first performed against the CGDB database for any input sequence or structure. Matches are subsequently realigned to the query using MUSCLE [77] for sequence realignment or “MAMMOTH” [78] for structure super-positioning. Residues in the query are then annotated as N (not in contact with the membrane), H (in contact with polar head group of the membrane lipids), or T (in contact with the lipid hydrophobic tails).

3.2. Tools for membrane MD simulation analysis

As MD simulations for membrane and membrane protein systems became widespread, many groups began developing tools to allow more efficient analysis of the MD trajectories. APL is an important indicator of the membrane phase and stability of the simulation. “GridMAT-MD” is a Perl program which can calculate the APL as well as thickness of a membrane [79]. For bilayer thickness calculation, the user can define a reference atom (such as the lipid phosphate, P atom) and the program first uses the upper leaflet as a reference and assign a paired lipid in the lower leaflet with the upper leaflet based on proximity in the x- and y-direction. The z-distance between the two points is calculated. The program then repeats the same step using the bottom leaflet as reference, and the two results are averaged, and written to a generic ASCII.dat file. Meanwhile, APL calculation of lipid-only systems can be as simple as taking the box size divided by the number of lipids in the upper or lower leaflet. Calculation of APL for membrane-protein systems are not as simple, and “GridMAT-MD” solves the problem by assigning protein atoms found within the lipid head groups to grid points then subtracts the total protein area from the size of the system. As of version 2.0, “GridMAT-MD” can now calculate the bilayer thickness and APL for multi “.pdb” or “.gro” files.

Mori and colleagues proposed a more sophisticated method for calculating the APL using Voronoi tessellation and Monte Carlo simulation [80]. Coordinates of center of mass for each lipid molecules and coordinates for protein atoms located between the maximum and minimum z-coordinates for the monolayer are projected onto the XY plane. Two-dimensional Voronoi analysis is subsequently performed for the lipids only. The APL for non-boundary lipids is the area of the Voronoi polygon where the lipid center of mass is located. The APL for boundary lipids can be determined by using a Monte Carlo integration method where the lipid region is probed by randomly shooting a pseudoparticle into the lipid Voronoi polygon. Thus, the APL for the boundary lipid is the product of the area of the Voronoi polygon, and the probability of the shot missing a protein atom. This method finds application in analysis of membrane-protein system and can differentiate between boundary and non-boundary lipids.

Not only that, analysis of other membrane properties became easier with the development of “Membrainy,” an intelligent membrane analysis tool that can provide the calculation of various membrane-specific properties for planar bilayer trajectories [81]. This include APL, order parameter, head group orientation, lipid mixing/demixing entropy, time evolution of the transmembrane voltage, 2D surface map generation, gel percentage, membrane thickness, detection of lipid flip-flop and annular shell lipid analysis. While the program has been primarily designed for use with Gromacs MD package, it is also compatible with pdb trajectory from other MD packages. Currently, it is implemented with CHARMM36, Berger/GROMOS87 and Martini v2.0 force fields, but is also expandable to include other force fields and trajectory formats. Output graphs can be readable by the Grace plotting software.

Other than that, MEMBPLUGIN [82] is another tool to study the MD trajectories of membrane-protein and complex membrane structures. This is a plugin in Visual MD package to measure biophysical properties in the simulated membranes.

4. Towards realistic bilayer simulations

Beyond the improvements in computational power, force field developments, and CG methodologies, more accurate representations of the membrane continued to evolve. The biological membrane is a complex entity composed of numerous lipid species such as phosphatidylcholines, phosphatidylethanolamines, and phosphatidylserine. Other molecules such as cholesterols, sphingomyelins, and cardiolipins also play a role in regulating membrane structure and function. In the outer membrane of many species of Gram-negative bacteria, the presence of lipopolysaccharide in the upperleaflet modulates the insertion, folding, and dynamics of outer membrane proteins within the membrane. Available tools for generate mixed membrane bilayers are including CHARMM-GUI Membrane Builder [73] and MemBuilder [83]. CHARMM-GUI Membrane Builder and MemBuilder supports a total of 32 and 18 different lipids types, respectively.

Straatsma and Soares first reported the simulation of the outer membrane protein OprF in an asymmetric outer membrane with a lipopolysaccharide and phospholipids, describing the saccharides component using GLYCAM parameters [84]. Holdbrook et al. performed simulations of the Haemphilus influenza Hia autotransporter domain in LPS and a realistic outer membrane inner leaflet which comprises 1-myristoyl 2-palmitoleoyl phosphatidylethanolamine (DMPE) lipid [85]. Comparison with simulations of the autotransporter in simpler, single species DMPC lipid model showed that the DMPC membrane accurately replicated the membrane thickness of the outer membrane and reproduced similar dynamics of the protein in asymmetric LPS/MPoPE membrane [85]. The realistic bilayer, however, revealed a patch of positive lysine and arginine residues on the extracellular mouth of Hia that interact regularly with phosphate and sugar groups of the LPS and are suggested to anchor Hia within the outer membrane [85].

5. Conclusion

The continuous update and improvement of atomistic force fields expanded the types of lipid molecules which could be simulated and increased the accuracy to better match experimental data. Depending on the level of detail of the simulation, UA force fields are an excellent alternative to balance between accuracy and speed. By using CG force fields approaches, sampling and size limitations may be tackled efficiently. Ultimately, advances in computational power and hardware have improved the timescale and system size where MD can be employed. In membrane simulations, the microsecond mark has been reached and simulations are slowly becoming routine work to complement experimental results. In addition, various web server tools and useful analysis programs have been developed to aid membrane simulation analysis. Further advances in lipid force fields will make it possible to characterize membrane structures in greater time and physical scale.

Acknowledgments

This work was funded by the Malaysia Ministry of Higher Education Fundamental Research Grant Scheme (FRGS) (203/CIPPM/6711439) and the Higher Institutions Center of Excellence (HICoE) Grant. S. W. Leong would also like to thank the Malaysia Ministry of Higher Education for MyBrain Science scholarship.