The authors thank Jinhui Li and Ricardo X. Ramirez for their simulation trajectories and discussions during the writing of this manuscript. O.C. was supported by the University at Buffalo Presidential Fellowship and National Institute of Health&#39;s Initiative for Maximizing Student Development Training Grant 1T32GM144920-01 awarded to Margarita L. Dubocovich (PI).

1. Building the system coordinates
<ol>
	<li>Navigate to CHARMM-GUI.org (C-GUI) using a web browser. On the top menu, navigate to Input Generator, then select Membrane Builder from the vertical options on the left side of the screen.</li>
	<li>To build a bilayer, select Bilayer Builder. 
	NOTE: First time users must activate their free account before building their first set of coordinates.</li>
	<li>Select Membrane Only System. Save the generated JOB ID to retrieve the system and continue from where you left off during the process if needed.
	<ol>
		<li>Visualize the systems during every step of the building process by clicking View Structure in the box located at the top of the page, or by downloading the resulting PDB file. Look out for missing components, mistakes in the selected input lipid species or patch size.</li>
	</ol>
	</li>
	<li>Select the components of the system.
	<ol>
		<li>Choose the Heterogeneous Lipid option, even if building a single-component bilayer; then select a Rectangular Box type.</li>
		<li>Select 45 water molecules per lipid for the hydration option; this is sufficient to ensure a fully hydrated bilayer.</li>
		<li>Set the length of XY to be based on number of lipid components. Then, select the number of lipids to be included for each lipid species determined ahead of the model. For the case study discussed in the next section, a membrane model with 600 lipids distributed symmetrically in two leaflets was built. To model the ER of eukaryotic cells a mixture of 336 DOPC, 132 DPPE, 60 CHOL, 72 POPI lipids was used for the PI model; and 330 DOPC, 126 DPPE, 54 CHOL, 66 POPI, and 24 DOPS lipids for the PI-PS model. 
		NOTE: C-GUI provides a library of lipid structures to choose from; click on the images next to the species name for its chemical structure.</li>
		<li>Type the desired number of molecules in the upper and lower leaflet in the two boxes next to the lipid name. For the case study, a symmetric membrane composition is desired - make sure there are no errors about unmatched number of lipids in the upper leaflet and lower leaflet. If asymmetry is desired, make sure the total number of lipids in each leaflet is correct. For details on building asymmetric bilayers, refer to work by Park et al.52,53.</li>
		<li>Go to the top of the lipid species list and click on Show the System Info&#160;button. Assemble components and complete the system.</li>
		<li>Select the option to include neutralizing ions using the distance-based algorithm for faster convergence54.</li>
		<li>Leave the default solution concentration of KCl at 0.15 mM. This is a&#160;typical salt concentration to render neutral the simulation box for membrane bilayers. 
		NOTE: If a different concentration is to be used, make sure to click the Calculate Solvent Composition button after editing it.</li>
	</ol>
	</li>
	<li>Select the simulation conditions and settings.
	<ol>
		<li>Select CHARMM36m as the FF option; it is commonly used for lipid and protein simulations, but the user may select other options discussed in the introduction.</li>
		<li>Select GROMACS as the MD engine to obtain sample input files in the corresponding format. 
		NOTE: GROMACS is recommended for new users because it has multiple online resources, tutorials, and forums for support. The user can select from multiple MD engines to explore options in terms of simulation performance and code syntax.</li>
		<li>Select the Constant Particle-Pressure-Temperature (NPT) ensemble, by far the most used dynamic ensemble in the simulation of lipid bilayers.</li>
		<li>Set the temperature and pressure in Kelvin and bars as 303 K and 1 bar, respectively. It is typical to set the temperature between 298 K and 310 K for the study of biological processes to ensure a bilayer in the liquid disorder state. 
		NOTE: The temperature depends on the conditions of the process to be simulated and can be altered as needed. Depending on the lipid species in the model, set the temperature to be above the transition temperature of pure lipid components prior to running the simulation.</li>
	</ol>
	</li>
	<li>Download the resulting files and transfer to the computer cluster.
	<ol>
		<li>Visualize the final system on a software of choice, such as VMD or PyMol, and inspect for proper setup. 
		&#8203;NOTE: It is good to check, for example, that there is sufficient water around the membrane such that lipids do not interact with image atoms during the simulation, and proper leaflet setup (a bilayer with no space or water in-between).</li>
	</ol>
	</li>
</ol>
2. Running MD simulations
<ol>
	<li>Upload and unzip the files from C-GUI on your computing cluster. Navigate to the Gromacs directory. Create a relaxation submission script.</li>
	<li>Follow the guidelines of the cluster for the format of a submission script.
	<ol>
		<li>Copy the commands listed until just above the # Production comment in the README file into the submission script. 
		NOTE: This default from C-GUI is a loop that runs a 6-step relaxation of the system. If a different and well-established protocol is desired, edit it to read the coordinates just built and downloaded from C-GUI.</li>
	</ol>
	</li>
	<li>Submit the relaxation script and confirm that all output files have been downloaded for all steps before moving to the production run. Upon completion, check for the following output files from GROMACS, generated during the 6 step run: *.log, *.tpr, *.gro, *.edr, *.trr / *.xtc</li>
	<li>Create a production run script.
	<ol>
		<li>Use one of the samples gmx grompp and gmx mdrun commands from any of the relaxation steps as a template.</li>
		<li>Before using the script, make sure to create a *.mdp file that contains similar simulation options as the provided step7_production.mdp file. 
		NOTE: The provided default options are standard for membrane simulations; ditances are listed in nm and time is given in picoseconds or number of steps (picoseconds / integration time step). Update the nsteps to run up to the desired simulation length (this is equal to dt * nsteps) and nst[x,v,f]out to update the data saving frequency in number of integration steps. For the case study, set the nsteps to 250,000,000 for a simulation length of 500ns (simulation time / integration step = 500,000 ps / 0.002ps), and nst[x,v,f]out to 50,000 to save data every 100 ps</li>
	</ol>
	</li>
	<li>Before running the actual simulation, run benchmark studies to determine the best use of resources.
	<ol>
		<li>Run the system for 1-2 ns using different number of computing nodes. 
		NOTE: The ER case study was submitted on the UB center for computational research (CCR) high-performance computing cluster55 for 2 ns, where the performance was tested for 1-10 nodes.</li>
		<li>Compare the performance in ns/day for each setting to determine the optimal resources for the run. It is common practice to select the number of nodes that result in 75%-80% of maximum performance.</li>
	</ol>
	</li>
	<li>Run the production run.
	<ol>
		<li>Run each system in triplicates to ensure reproducibility and conduct statistical analysis on the data.</li>
		<li>Extend the trajectory based on the benchmarks if the allowed queue time for submission on the computing cluster runs out. Use the gmx convert-tpr, then gmx mdrun commands to continue trajectory collection. 
		&#8203;NOTE: Options are described in the GROMACS documentation online (https://manual.gromacs.org/).</li>
		<li>For a membrane-only system, examine if the system has reached equilibrium by calculating the area-per-lipid over time. If it has not, extend the simulation trajectory.</li>
	</ol>
	</li>
</ol>
3. Analyzing the trajectory
<ol>
	<li>Visualize the system before running analysis to determine molecules of interest and portion of the trajectory that is intended for characterization.</li>
	<li>Compress raw trajectory files (*.trr) by changing file format to *.xtc and/or skipping frames to reduce file size and facilitate more efficient transfer to the local station for visualization and analysis. 
	NOTE: For large membrane systems, one may opt to strip water from the trajectory to further reduce the file size. This can be done with index files on GROMACS, TCL scripts on VMD, or Python libraries such as MDAnalysis and MDTraj.</li>
	<li>Perform chosen analyses during the equilibrated portion of the trajectory, as determined from the area-per-lipid timeseries. 
	NOTE: Refer to discussion for more details on typical membrane analyses and how to run them.</li>
</ol>

Simulating and Analyzing Lipid Bilayers Using Molecular Dynamics

The authors have no competing interests to disclose.

Experimental techniques can visualize biomolecules at high resolution using cryo-electron microscopy (cryo-EM)58, fluorescence techniques, and atomic force microscopy (AFM)59. However, it is challenging to capture the interplay and dynamics of molecular interactions that underlie biological pathways, disease pathogenesis, and therapeutic delivery at the atomic or amino acid level. Here, capabilities of MD simulations to study lipid membranes and the main steps to design, bu...

Lipids are major constituents of membranes, which provide boundaries for cells and enable intracellular compartmentalization1,2,3. Lipids are amphiphilic, with a polar head group and two hydrophobic fatty acid tails; these self-assemble into a bilayer to minimize contact of the hydrophobic chains with water3,4. Various combinations of hydrophilic head groups and hydrophobic tails result in different classes of lipids in biological membranes, such as glycerophospholipids, sphingolipids, and sterols (Figure 1)1,5,6. Glycerophospholipids are primary building blocks of eukaryotic cell membranes composed of glycerophosphate, long-chain fatty acids, and head groups of low-molecular weight7. Lipid nomenclature is based on differences in head groups; examples include phosphatidyl-choline (PC), phosphatidyl-ethanolamine (PE), phosphatidyl-serine (PS), phosphatidyl-glycerol (PG), phosphatidyl-inositol (PI), or the unmodified phosphatidic acid (PA)5,6. As for hydrophobic tails, the length and degree of saturation vary, along with the backbone structure. The possible combinations are numerous, resulting in thousands of lipid species in mammalian cells6. Changes in membrane lipid composition lead to different mechanical and structural membrane properties that impact the activity of both integral membrane proteins and peripheral proteins2,6.
<img alt="Figure 1" class="xfigimg" src="/files/ftp_upload/65712/65712fig01.jpg" /> 
Figure 1. Representative lipid structures. Fatty acid tails are shown in blue boxes, common lipid head groups in orange, and sample backbones in purple. <a href="https://www.jove.com/files/ftp_upload/65712/65712fig01large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
Lipids are active players in cellular processes, protein activation in signaling cascades, and healthy cell homeostasis8,9. Altered lipid dynamics are the result of infection or can be markers of pathogenesis of disease10,11,12,13,14,15. As barriers for the cell, the study of membrane lipids and their role in permeation of small molecules is of relevance for drug delivery systems and membrane disruption mechanisms16,17. Chemical diversity and different ratios of lipid species across organelles, tissues, and organisms give rise to complex membrane dynamics2. It is therefore important to retain these characteristics in modelling studies of lipid bilayers, especially when the goal of a study is to examine interactions of other biomolecules with the membrane. The lipid species to consider in a model depends on the organism and cellular compartment of interest. For example, PG lipids are important for electron transfer in photosynthetic bateria18, while phosphorylated inositol lipids (PIPs) are major players in plasma membrane (PM) dynamics and signaling cascades in mammalian cells19,20. Inside the cell, the PM, endoplasmic reticulum (ER), Golgi, and mitochondrial membranes contain unique lipid abundances that influence their function. For example, the ER is the hub for lipid biogenesis and transports cholesterol out to the PM and Golgi; it contains a high lipid diversity with an abundance of PC and PE, but low sterol content, which promotes membrane fluidity21,22,23,24. In contrast, the PM incorporates hundreds and even thousands of lipid species depending on the organism25, it contains high levels of sphingolipids and cholesterol that give it a characteristic rigidity compared to other membranes in the cell24. Leaflet asymmetry should be considered for membranes like the PM, which has an outer leaflet rich in sphingomyelin, PC, and cholesterol, and an inner leaflet rich in PE, PI and PS that are important for signaling cascades24. Finally, lipid diversity also prompts the formation of micro-domains that differ in packing and internal order, known as lipid rafts24,26; these exhibit lateral asymmetry, are hypothesized to play important roles in cellular signaling26, and are hard to study due to their transient nature.
Experimental techniques such as fluoroscopy, spectroscopy, and model membrane systems like giant unilamellar vesicles (GUVs) have been used to investigate interactions of biomolecules with membranes. However, the complex and dynamic nature of the components involved is difficult to capture with experimental methods alone. For example, there are limitations on the imaging of transmembrane domains of proteins, the complexity of membranes used in such studies, and the identification of intermediate or transient states during the process of interest27,28,29. Since the advent of molecular simulation of lipid monolayers and bilayers in the 1980s29, lipid-protein systems and their interactions can now be quantified at the molecular level. Molecular dynamics (MD) simulation is a common computational technique that predicts the movement of particles based on their intermolecular forces. An additive interaction potential describes the bonded and non-bonded interactions between particles of the system30. The set of parameters used to model these interactions is termed the simulation forcefield (FF). These parameters are obtained from ab initio computations, semi-empirical, and quantum mechanical calculations, and optimized to reproduced data from X-ray and electron diffraction experiments, NMR, infrared, Raman and neutron spectroscopy, amongst other methods31.
MD simulations can be used to study systems at various levels of resolution32,33,34. Systems that aim to characterize specific biomolecular interactions, hydrogen bonds, and other high-resolution details are studied with all-atom (AA) simulations. In contrast, coarse-grained (CG) simulations lump atoms into larger functional groups to reduce computational cost and examine larger scale dynamics33. Situated in between these two are united-atom (UA) simulations, where hydrogen atoms are combined with their respective heavy atoms to accelerate the computation33,35. MD simulations are a powerful tool for exploration of the dynamics of lipid membranes and their interactions with other molecules and can serve to provide molecular level mechanisms for processes of interest at the membrane interface. Additionally, MD simulations can serve to narrow down experimental targets and predict macromolecular properties of a given system based on microscopic interactions.
In brief, given a set of initial coordinates, velocities, and a set of conditions like constant temperature and pressure, positions and velocities of each particle are calculated through numerical integration of the interaction potential and Newton&#39;s Law of motion. This is repeated iteratively, thereby generating a simulation trajectory30. These calculations are performed with an MD engine; among several open-source packages, GROMACS36 is one of the most commonly used engines and the one we describe here. It also includes tools for analysis and constructing initial coordinates of systems to be simulated37. Other MD engines include NAMD38; CHARMM39, and AMBER40, which the user can select at their own discretion based on computational performance of a given system. It is critical to visualize the trajectories during the simulation as well as for analysis and interpretation of the results. A variety of tools are available; here we discuss visual molecular dynamics (VMD) that offers a wide range of features, including three dimensional (3-D) visualization with expansive drawing and coloring methods, volumetric data visualization, building, preparing, and analyzing trajectories of MD simulation systems, and trajectory-movie making with no limits on system size, if the memory is available41,42,43.
The accuracy of predicted dynamics between system components is directly influenced by the FF chosen for the propagation of the trajectory. Empirical FF parametrization efforts are pursued by few research groups. The most established and common FF for MD include CHARMM39, AMBER40, Martini44, OPLS45, and SIRAH46. The all-atom additive CHARMM36 (C36) force field47 is widely used for AA MD of membrane systems as it accurately reproduces experimental structural data. It was originally developed by the CHARMM community, and it is compatible with multiple MD engines like GROMACS and NAMD. Despite improvements across common FFs, there is a continuous effort to improve the parameter sets to allow predictions that closely reproduce experimental observables, driven by interests in particular systems of study48,49.
A challenge when simulating lipid membranes is determining the length of simulation trajectory. This is largely dependent on the metrics to be analyzed and the process that one aims to characterize. Typically, complex lipid mixtures require longer time to reach equilibrium, as more species must have enough time to diffuse on the membrane plane and reach a stable lateral organization. A simulation is said to be at equilibrium when the property of interest has reached a plateau and fluctuates about a constant value. It is common practice to obtain at least 100-200 ns of equilibrated trajectory to carry out appropriate statistical analysis on the properties and interactions of interest. It is common to run membrane-only simulations between 200-500 ns, depending on the complexity of the lipid mixture and research question. Protein-lipid interactions typically require longer simulation times, between 500-2000 ns. Some approaches to accelerate sampling and observable dynamics with membrane systems are: (i) the highly mobile membrane mimetic (HMMM) model, which replaces end carbons of lipids in the membrane with organic solvent to accelerate sampling50; and (ii) hydrogen mass repartitioning (HMR), which combines a fraction of the masses of heavy atoms within a system with those of hydrogen atoms to allow the use of a larger simulation timestep51.
The following protocol discusses a beginner-friendly approach to build, run, and analyze realistic membrane models using AA MD. Given the nature of MD simulations, multiple trajectories must be run to account for reproducibility and proper statistical analysis of the results. It is current practice to run at minimum three replicas per system of interest. Once the lipid species have been selected for the organism and process of interest, basic steps to build, run, and analyze a simulation trajectory of a membrane-only system are outlined and summarized in Figure 2.
<img alt="Figure 2" class="xfigimg" src="/files/ftp_upload/65712/65712fig02.jpg" /> 
Figure 2. Schematic to run MD simulations. Orange boxes correspond to the three main steps described in the protocol. Underneath is the workflow of the simulation process. During system setup, the system containing the initial coordinates of a solvated membrane system is built with a system input generator like CHARMM-GUI Membrane Builder. After transferring the input files to a high-performance computing cluster, the simulation trajectory is propagated using an MD engine, such as GROMACS. Trajectory analysis can be done on the computer cluster or a local working station along with visualization. Analysis is then carried out, using either packages with built-in analysis code such as GROMACS and VMD, or by using Bash scripts or various Python libraries. <a href="https://www.jove.com/files/ftp_upload/65712/65712fig02large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>

<table border="1" fo:keep-together.within-page="1" fo:keep-with-next.within-page="always">
	<tbody>
		<tr>
			<td>Anaconda3</td>
			<td>Anaconda Inc (Python &#38; related libraries)</td>
			<td>N/A</td>
			<td></td>
		</tr>
		<tr>
			<td>CHARMM-GUI.org</td>
			<td>Im lab, Lehigh University</td>
			<td>N/A</td>
			<td></td>
		</tr>
		<tr>
			<td>GROMACS</td>
			<td>GROMACS development team</td>
			<td>N/A</td>
			<td></td>
		</tr>
		<tr>
			<td>Linux HPC Cluster</td>
			<td>UB CCR</td>
			<td>N/A</td>
			<td></td>
		</tr>
		<tr>
			<td>MATLAB</td>
			<td>MathWorks</td>
			<td>N/A</td>
			<td></td>
		</tr>
		<tr>
			<td>VMD</td>
			<td>Theoretical and Computational Biophysics Group</td>
			<td>N/A</td>
			<td></td>
		</tr>
	</tbody>
</table>

realistic membrane modeling using complex lipid mixtures in simulation studies

Lipids are structural building blocks of cell membranes; lipid species vary across cell organelles and across organisms. This variety results in different mechanical and structural properties in the membrane that directly impact the molecules and processes that occur at this interface. Lipid composition is dynamic and can serve to modulate cell signaling processes. Computational approaches are increasingly used to predict interactions between biomolecules and provide molecular insights to experimental observables. Molecular dynamics (MD) is a technique based on statistical mechanics that predicts the movement of atoms based on the forces that act on them. MD simulations can be used to characterize the interaction of biomolecules. Here, we briefly introduce the technique, outline practical steps for beginners who are interested in simulating lipid bilayers, demonstrate the protocol with beginner-friendly software, and discuss alternatives, challenges, and important considerations of the process. Particularly, we emphasize the relevance of using complex lipid mixtures to model a cell membrane of interest to capture the appropriate hydrophobic and mechanical environments in simulation. We also discuss some examples where membrane composition and properties modulate the interactions of bilayers with other biomolecules.

Author Spotlight: Advancing Cell Membrane Biophysics - Exploring Interactions and Challenges Through Experimental and Computational Approaches 

Realistic Membrane Modeling Using Complex Lipid Mixtures in Simulation Studies

To illustrate the use of the protocol and the results that can be obtained, a comparison study for membrane models for the endoplasmic reticulum (ER) is discussed. The two models in this study were (i) the PI&#160;model, which contains the top four lipid species found in the ER, and (ii) the PI-PS model, which added the anionic phosphatidylserine (PS) lipid species. These models were later used in a study of a viral protein and how it interacts with the membrane, the interest on PS has been cited as important for the per...

Watch this Scientific Journal Video about Advancing Cell Membrane Biophysics - Exploring Interactions and Challenges Through Experimental and Computational Approaches at JoVE.com

Membrane lipid diversity in structure and composition is an important contributor to cellular processes and can be a marker of disease. Molecular dynamics simulations allow us to study membranes and their interactions with biomolecules at atomistic resolution. Here, we provide a protocol to build, run, and analyze complex membrane systems.

Lipids are structural building blocks of cell membranes; lipid species vary across cell organelles and across organisms. This variety results in different ...

realistic-membrane-modeling-using-complex-lipid-mixtures-simulation

Research

JoVE Journal

Biology

1.6K Views.  State University of New York at Buffalo. Membrane lipid diversity in structure and composition is an important contributor to cellular processes and can be a marker of disease. Molecular dynamics simulations allow us to study membranes and their interactions with biomolecules at atomistic resolution. Here, we provide a protocol to build, run, and analyze complex membrane systems.