Martini tutorials: lipids with lipidome
- Details
- Created: Tuesday, 20 May 2014 08:34
- Last Updated: Thursday, 11 October 2018 13:01
Martini tutorials: lipids with the lipidome
Summary
Advanced: refine lipid parameters
Tools and scripts used in this tutorial
Lipid bilayers
The Martini coarse-grained (CG) model was originally developed for lipids[1,2]. The underlying philosophy of Martini is to build an extendable CG model based on simple modular building blocks and to use only few parameters and standard interaction potentials to maximise applicability and transferability. Martini uses an approximate 4:1 atomistic (heavy atoms) to CG bead mapping and in version 2.0[2] 18 bead types were defined to represent chemical characteristics of the beads. The CG beads have a fixed size and interact using an interaction map with 10 different interaction strengths. Due to the modularity of Martini, a large set of different lipid types has been parameterized. Their parameters are available in the Martini Lipidome, see the Martini lipid topology section. In this tutorial, you will learn how to set up lipid-water system simulations with the lipidome, with a focus on bilayers. You will also study a number of standard bilayer properties.
These tutorials assume a basic understanding of the Linux operating system and the gromacs molecular dynamics (MD) simulation package. An excellent gromacs tutorial is available at: cgmartini.nl/~mdcourse/.
The aim of the tutorial is to create and study properties of CG Martini models of lipid bilayer membranes. First, we will attempt to spontaneously self-assemble a DPPC bilayer and check various properties that characterize lipid bilayers, such as the area per lipid, bilayer thickness, order parameters, and diffusion. Next, we will change the nature of the lipid head groups and tails, and study how that influences the properties. Finally, we will move on to creating more complex multicomponent bilayers.
Download all the files of this tutorial: bilayer_lipidome_tutorial.tgz Unpack the lipidome-tutorial.tgz archive (it expands to a directory called bilayer-lipidome-tutorial), and enter the bilayer-lipidome-tutorial directory:
$ tar -xzvf bilayer_lipidome_tutorial.tgz
$ cd bilayer-lipidome-tutorial
Bilayer self-assembly
We will begin with self-assembling a dipalmitoyl-phosphatidylcholine (DPPC) bilayer from a random configuration of lipids and water in the simulation box. Enter the spontaneous-assembly subdirectory. The first step is to create a simulation box containing a random configuration of 128 DPPC lipids. This can be done by starting from a file containing a single DPPC molecule. A file with coordinates for a single DPPC molecule is available for you as dppc_single.gro. The gromacs tool genbox can take this single-molecule conformation and attempt to place it in a simulation box multiple times at a random position and random orientation, each time checking that there are no overlaps between the consecutively placed molecules.
$ genbox -ci dppc_single.gro -nmol 128 -box 7.5 7.5 7.5 -try 500 -o 128_noW.gro
Perform a short energy minimization of the system containing only the lipids; the only reason for doing this, is getting rid of high forces between beads that may have been placed quite close to each other. The settings file minimization.mdp is provided for you, but you will need the topology for the DPPC lipid. Download it from the Martini lipidome page, see the Martini lipid topology section, place it in your current working directory, and check that the file name corresponds to the #include statement in the topology file dppc.top.
$ grompp -f minimization.mdp -c 128_noW.gro -p dppc.top -o dppc-min-init.tpr
$ mdrun -deffnm dppc-min-init -v -c 128_minimized.gro
Using the gromacs tool genbox again, add 6 CG waters (note that this corresponds to 24 all-atom waters) per lipid.
$ genbox -cp 128_minimized.gro -cs water.gro -o waterbox.gro -maxsol 768 -vdwd 0.21
The value of the flag -vdwd (default van der Waals radii) has to be increased from its default atomistic length (0.105 nm) to a value reflecting the size of Martini CG beads. The water beads are taken from the file water.gro provided for you. A new file, waterbox.gro is produced containing the 128 lipids and added water beads.
Now perform an energy minimization again. You will need to adapt dppc.top to reflect the water (W) beads added to the system. Look at the output of genbox to see how many W beads were added; the maximum number is 768 (specified on the last command as the -maxsol option). Remember that in gromacs files a semicolon (;) preceeds a comment line and is not interpreted as input!
$ [gedit/vi/other editor] dppc.top
$ grompp -f minimization.mdp -c waterbox.gro -p dppc.top -o dppc-min-solvent.tpr
$ mdrun -deffnm dppc-min-solvent -v -c minimized.gro
Now, you are ready to run the self-assembly MD simulation. About 30 ns should suffice.
$ grompp -f martini_md.mdp -c minimized.gro -p dppc.top -o dppc-md.tpr
$ mdrun -deffnm dppc-md -v
This might take approximately 45 minutes on a single CPU but by default mdrun will use all available CPUs on your machine. You can tune the numbers of parallel threads mdrun uses with the -nt parameter. You may want to check the progress of the simulation to see whether the bilayer has already formed before the end of the simulation. You may do this most easily by starting the inbuilt gromacs viewer (you will need to open a new terminal and make sure you are in the directory where this simulation is running):
$ ngmx -f dppc-md.xtc -s dppc-md.tpr
or, alternatively, use VMD or pymol. Both VMD and pymol suffer from the fact that Martini bonds are usually not drawn because they are much longer than standard bonds and the default visualization is not very informative. For visualization with pymol this is solved most easily by converting the trajectory to .pdb format, explicitly writing out all bonds. The disadvantage is that very large files are produced in this conversion!
$ trjconv -s dppc-md.tpr -f dppc-md.xtc -o dppc-md.pdb -pbc whole -conect
$ pymol dppc-md.pdb
For VMD, a plugin script cg_bonds.tcl was written that takes the Gromacs topology file and adds the Martini bonds defined in the topology. In general, a useful preprocessing step for VMD visualization is to avoid molecules being split over the periodic boundary, because if they are, very long bonds will be drawn between them. A script do_vmd.sh has been prepared for you for visualization using VMD. To learn more about visualizing Martini systems using VMD, consult the VMD tutorial.
$ ./do_vmd.sh
The inital and final snapshots should look similar to Fig. 1, at least if the self assembly resulted in a bilayer in the alloted time, which is not guaranteed.
In the meantime, have a close look at the Martini parameter files. The available bead types and their interactions are defined in the file martini_v2.1.itp. This file contains the definitions of version 2.1, and also includes the definition of the water bead, right at the bottom of the file. The lipid DPPC is defined in the file martini_v2.0_DPPC_01.itp , which you downloaded from the lipidome page. Get out the 2007 Martini paper[2] and check the definition of this molecule. Understanding how these files work, will help you work with the Martini model and define new molecules or refine existing models.
Figure 1 | DPPC bilayer formation. A) 64 randomly localized DPPC lipid molecules with water. B) After a 30 ns long simulation the DPPC lipids have aggregated together and formed a bilayer.
When the simulation has finished, check whether you got a bilayer. If yes, check if the formed membrane is normal to the z-axis (i.e., that the membrane is situated in the xy-plane). Have a look at the self-assembly process: can you recognize intermediate states, such as micelles, and do you see metastable structures such as a water pore and/or a lipid bridge?
Bilayer equilibrium run and analysis
If your bilayer was formed in a plane other than the xy-plane, rotate the system so that it will (with editconf). In case you did not get a bilayer at all, you can continue the tutorial with the pre-formed one from dspc-bilayer.gro. From this point on, we will refer to the bilayer as a DSPC bilayer. In Martini, the 4-to-1 mapping limits the resolution and the same model represents in this case a C-16 tail (palmitoyl) and a C-18 tail (stearoyl, hence DSPC). The reason for referring to the lipid as DSPC is that later on, we will compare the simulated properties to experiment and we will compare different lipids all with the same number of beads, changing either the head group from PC to PE or changing the tail from stearoyl (C-18, fully saturated) to oleoyl (C-18, single cis double bond at C-9).
The spontaneous assembly simulation was done using isotropic pressure coupling. The bilayer may have formed but is probably under tension because of the isotropic pressure coupling. Therefore, we first need to run a simulation in which the area of the bilayer can reach a proper equilibrium value. This requires that we use independent pressure coupling in the plane and perpendicular to the plane. Set up a simulation for another 30 ns at zero surface tension (switch to semiisotropic pressure coupling; if the pressure is the same in the plane and perpendicular to the plane, the bilayer will be at zero surface tension) at a higher temperature of 341 K. This is above the experimental gel to liquid-crystalline transition temperature of DSPC. You will find how to change pressure coupling and temperature in the gromacs manual: manual.gromacs.org/current/online/mdp_opt.html Again, you may not want to wait for this simulation to finish. In that case, you may skip ahead and use a trajectory provided in the subdirectory spontaneous-assembly/equilibration/.
From the latter simulation, we can calculate properties such as:
- area per lipid
- bilayer thickness
- P2 order parameters of the bonds
- lateral diffusion of the lipids
In general, for the analysis, you might want to discard the first few ns of your simulation as equilibration time. The area per lipid as a function of time should give you an indication of which part of the simulation you should definitely not include in the analysis.
Area per lipid
To get the (projected) area per lipid, you can simply divide the area of your simulation box by half the number of lipids in your system. The box-sizes can be obtained by running the gromacs tool g_energy (ask for Box-X and Box-Y). Then either use awk or xmgrace to calculate the area per lipid as a function of time. The gromacs tool g_analyze provides a convenient way to calculate the average and error estimate in any series of measurements (use the option -ee for the error estimate, and make sure you feed g_analyze a file with a .xvg extension). (Note that this calculation might not be strictly OK, because the self-assembled bilayer might be slightly asymmetric in terms of number of lipids per monolayer, i.e., the areas per lipid are different for the two monolayers. However, to a first approximation, we will ignore this here.)
Bilayer thickness
To get the bilayer thickness, use g_density. You can get the density for a number of different functional groups in the lipid (e.g., phosphate and ammonium headgroup beads, carbon tail beads, etc) by feeding an appropriate index-file to g_density (make one with make_ndx; you can select, e.g., the phosphate beads by typing a P*). The bilayer thickness you can obtain from the distance between the headgroup peaks in the density profile.
A more appropriate way to compare to experiment is to calculate the electron density profile. The g_density tool also provides this option. However, you need to supply the program with a data-file containing the number of electrons associated with each bead (option -ei electrons.dat). The format is described in the gromacs manual.
Compare your results to those from small-angle neutron scattering experiments[3]:
- thickness = 4.98 ± 0.15 nm
- area per lipid = 0.65 ± 0.05 nm2
Lateral diffusion
Before calculating the lateral diffusion, remove jumps over the box boundaries (trjconv -pbc nojump). Then, calculate the lateral diffusion using g_msd. Take care to remove the overall center of mass motion (-rmcomm*), and to fit the line only to the linear regime of the mean-square-displacement curve (-beginfit and -endfit options of g_msd). To get the lateral diffusion, choose the -lateral z option.
To compare the diffusion coefficient obtained from a Martini simulation to a measured one, a conversion factor of about 4 has to be applied to account for the faster diffusion at the CG level due to the smoothened free energy landscape (note, the conversion factor can vary significantly depending on the molecule in question).
*In gromacs v3, this option is not available. An additional trjconv -center will do the trick.
Order parameters
Now, we will calculate the (second-rank) order parameter, which is defined as:
P2 = 1/2 (3 cos2<θ> − 1),
where θ is the angle between the bond and the bilayer normal. P2 = 1 means perfect alignment with the bilayer normal, P2 = −0.5 anti-alignment, and P2 = 0 random orientation.
A script to calculated these order parameters can be downloaded here (there is also a version located in the directory scripts/). Copy the .xtc and .tpr files to a analysis subdirectory (the script expects them to be named traj.xtc and topol.tpr, respectively; you may want to use the 30 ns simulations already present in the subdirectory spontaneous-assembly/equilibration/). The script do-order.py will calculate P2 for you. As it explains when you invoke it, it needs a number of arguments. The command:
$ ./do-order.py my.xtc my.tpr 15000 30000 20 0 0 1 128 DPPC
will for example read 15 to 30 ns from the trajectory my.xtc. The simulation has 128 DPPC lipids and the output is the P2, calculated relative to the Z-axis, and average over results from every 20th equally-spaced snapshot in the trajectory. (Remember, DSPC and DPPC share the same Martini representation!)
Advanced: additional analysis
Different scientific questions require different methods of analysis and no set of tools can span them all. There are various tools available in the gromacs package, see the gromacs manual: manual.gromacs.org/current/online.html. Most simulation groups, therefore, develop sets of customized scripts and programs many of which they make freely available, such as the Morphological Image Analysis and the 3D pressure field tools available here. Additionally a number of packages are available to for assistance with analysis and the development of customized tools, such as the python MDAnalysis package pypi.python.org/pypi/MDAnalysis/.
Changing lipid type
Lipids can be thought of as modular molecules. In this section, we investigate the effect of changes in the lipid tails and in the headgroups on the properties of lipid bilayers using the Martini model. We will i) introduce a double bond in the lipid tails, and ii) change the lipid headgroups from phosphatidylcholine (PC) to phosphatidylethanolamine (PE).
Unsaturated tails
To introduce a double bond in the tail, we will replace the DSPC lipids by DOPC. Set up a directory for the DOPC system. Enter that directory. Download the topology for DOPC from the Martini lipidome web page. Compare the Martini topologies of these two lipids. You will see that the number of beads for these lipids is the same. You can therefore set up a DOPC bilayer quite simply from your DSPC result. Copy over the .top and .mdp file from the equilibration DSPC run. Copy the final frame of the DSPC run to serve as the starting frame of the DOPC run. Replace DPPC by DOPC in your .top and .mdp files, and grompp will do the rest for you (you can ignore the ”atom name does not match” warnings of grompp). Execute the simulation, or - if you are impatient - use the trajectory provided for you as spontaneous-assembly/equilibration/dopc/dopc-ext.xtc. NOTE: because the lipids DSPC and DOPC have exactly the same number of beads, you can use the .tpr file from DSPC to analyze the dopc-ext.xtc file.
Changing the headgroups
Similarly, starting again from the final snapshot of the DSPC simulation, you can easily change the head groups from PC to PE. Download the topology for DPPE from the Martini lipidome web page. Also for this system, run a 30 ns MD simulation, or - if you are impatient - use the trajectory provided for you as spontaneous-assembly/equilibration/dspe/dspe-ext.xtc. NOTE: because the lipids DSPC and DOPC have exactly the same number of beads, you can use the .tpr file from DSPC to analyze the dopc-ext.xtc file.
Compare the above properties (area per lipid, thickness, order parameter, diffusion) between the three bilayers (DSPC, DOPC, DSPE). Do the observed changes match your expectations? Why/why not? Discuss.
Complex bilayer mixtures
When setting up larger and more complex bilayers it can be more convenient, or necessary, to start with a bilayer close to equilibrium rather then trying to get the bilayer to self-assemble. This can be done by concatenating (e.g. using editconf) and/or altering previously formed bilayers or by using a bilayer formation program such as insane.py, available in the directory scripts/ or downloadable here. INSANE (INSert membrANE) is a CG building tool that generates membranes by distributing lipids over a grid. Lipid structures are derived from simple template definitions that can either be added to the program or provided on-the-fly from combinations of headgroup, linker and lipid tails specified on the command-line. The program uses two grids, one for the inner and one for the outer leaflet, and distributes the lipids randomly over the grid cells with numbers according to the relative quantities specified. This allows for building asymmetric bilayers with specific lipid compositions. The program also allows for adding solvent and ions, using a similar grid protocol it distributes them over a 3D grid. Additional information about the functionality of INSANE can be found by running insane.py -h.
In the directory complex-bilayers/ we will create a fully hydrated 4:3:3 mole ratio bilayer of DPPC:DIPC:CHOL with physiological ion concentration using insane.py, this is a similar mixture as was used in Risselada and Marrink[4] to show raft formation. The lipid DIPC is representative for tails containing two cis-unsaturated bonds and of length 16-18 carbons. Start by running INSANE (downloaded or copied from the scripts/ directory):
$ python insane.py -l DPPC:4 -l DIPC:3 -l CHOL:3 -salt 0.15 -x 15 -y 10 -z 9 -d 0 -pbc cubic -sol W -o dppc-dipc-chol-insane.gro
This will generate an initial configuration for the system dppc-dipc-chol-insane.gro (Fig. 2A). Download the relevant lipid topologies from the Martin lipdome and provide this information together with the system composition to the .top file. (Check that the file provided for you is correct.)
Then, energy minimize the structure and gently equilibrate the system. Note, as this simulation contains multiple components you will have to make an index file (using make_ndx) and group all the lipids together and all the solvent together to fit the specified groups in the .mdp files; the .mdp files provided expect the names Lipids and Solvent. As all the bilayer lipids and solvent were placed on a grid (Fig. 2A), even after minimization they can still be in an energetically unfavorable state. Due to the large forces involved it may be necessary to run a few equilibrium simulations using a short time step (1-10 fs) before running production simulations with Martini lipid time steps (30-40 fs). The initial grid order imposed by insane.py should relax in a few ns (Fig. 2B), we recommend simulating for 5-30 ns using the Berendsen pressure coupling algorithm, to relax the membrane area, before switching to Parrinello-Rahman for the production run. This mixture will phase separate at a temperature of 295 K but that can take multiple microseconds (Fig. 2C). A sample run is provided for you in the directory bilayer-lipidome-tutorial/complex-bilayers/raftmix.xtc. The INSANE script can also be used to set up a complex (or simple!) bilayer system including membrane protein, see the Martini protein tutorial.
Figure 2 | DPPC-DIPC-CHOL bilayer formation. A) Using insane.py a 4:3:3 DPPC-DIPC-CHOL bilayer was setup. B) After a 20 ns long simulation the grid structure is gone and the bilayer area has equilibrated. C) This lipid mixture phase separates into Ld and Lo domains at 295 K, this phase separation can take multiple microseconds, showing a top view of the simulation after 6 microseconds. In all panels DPPC is in red, DIPC in blue and cholesterol in orange.
Advanced: refine lipid parameters
This section uses the backward.py tool for fine-grained-to-coarse-grained (FG-to-CG) transformation, which is treated in more detail elsewhere. A related, but different strategy not using the backward.py tool is used in the tutorial on polymers, advantageous there because the PEG molecule is simpler to map from FG to CG representation than lipids.
In this part, we will outline how atomistic simulations can be used to obtain improved Martini parameters for selected degrees of freedom. You may then take this strategy to whatever lengths you deem fit. Remember, Martini in general does not aim to be specific for one particular molecule, but recognizes that the four-to-one mapping results in Martini molecules being representative of several closely related molecules. However, you may want to develop a specific model or compare how well Martini covers closely related molecules.
Here, we study diC18:2-PC, a PC lipid with a doubly unsaturated tail, in the Martini lipidome denoted as DIPC. Specifically, the atomistic (FG) model is for dilineoylphosphatidylcholine. The lineoyl tail is fully denoted as C18:2c9,12. This notation shows that the tail is 18 carbons long (counting starts at the ester C-atom, bound to the two O-atoms), contains two cis bonds, starting at carbon 9 and carbon 12, respectively (i.e. the cis-bonds are between carbons 9-10 and 12-13). We will show how to compare distributions of selected angles between beads for a mapped FG model to those from a Martini simulation. You can then use the distributions to guess and refine selected Martini parameters such that the distributions match those from the FG simulation as close as possible. The files needed in this section are located in the bilayer-lipidome-tutorial/refine/ directory. The FG trajectory is in fg.xtc. This contains only the lipids from a simulation of a lipid bilayer consisting of 128 lipids (64 per monolayer) in water.
The task can be divided into the following five steps:
- Transform the all-atom (fine-grained, FG) trajectory into its coarse-grained counterpart, based on a pre-defined mapping scheme.
- Calculate selected angle distributions. These will serve as the reference (“gold standard”) later.
- Do a coarse-grained simulation of the system.
- Calculate the selected angle distibutions, and compare to the reference.
- Revise coarse-grained parameters, repeat steps 3 and 4 until you are satisfied with the result.
For the transformation from FG to CG, we will use the tool backward.py, which is described in more detail elsewhere. You do not need to download this tool; it is provided in the refine/ directory. The refine/Mapping directory contains a file dipc.gromos.map that defines the mapping from the FG to the CG model. Have a look at it. You will see how each (united) atom of the FG lipid contributes to a CG bead. The order of the CG beads in the CG model is provided under the [ martini ] directive. Beads will appear in the output file in this order. It is important that this order corresponds to your proposed Martini model. Check the correspondence by inspecting the file refine/martini_v2.0_DIPC.itp. This contains the proposed Martini model. If you wish to alter the mapping scheme, you can do so by assigning the atoms differently to the beads. In general, we aim to adapt the mapping such that we get predominantly unimodal distributions for bonds and angles; however, we also want to stay within the Martini scheme of 4-to-1 mapping and mapping to recognizable chemical groups. As long as you can rationalize your choices, Martini is quite flexible.
1. Transformation of the FG trajectory to a mapped CG trajectory
Backward.py can perform the mapping from the FG to the CG model for .pdb and .gro files, but not (yet) for entire trajectories. Thus, the FG trajectory is first written out in single frames, and each frame is mapped from FG to CG using backward.py. When all CG frames are done, they are concatenated into a single CG trajectory - the time stamp is lost in this process! The mapping procedure is done for you in the bash script do-mapping.sh. Please inspect it before you execute it. Before you do the mapping, you need to have a .tpr file for the FG system.
$ grompp -p system.top -c lipid.gro -f min.mdp -o fg.tpr -maxwarn 10
$ ./do-mapping.sh
You will see that two directories have been created: FRAMES contains all the FG frames dumped from the trajectory; CGFRAMES contains all the CG frames converted from the FG frames and the mapped trajectory (CGFRAMES/cg-mapped.xtc).
2. Calculation of the selected angle distributions
Angle distributions can be calculated using the gromacs tool g_angle. This requires an index file containing multples of three indices. The three indices are the atoms/beads that define an angle. In this case, you have 128 lipids that each contain one or more angles you are interested in. You want to pool the data from equivalent angles and collected over the entire trajectory. This is best done by calculating the angle of interest for each lipid as a function of time, collect all equivalent angle data in one file and generate a distribution from that data. The bash script angledistr.sh does that for you. However, before you start, you need to provide a list of angles you are interested in. The script expects a file called anglelist containing one or more lines consisting of three numbers. The numbers are the indices of the angle(s) of interest. The numbers refer to a single lipid. The script will convert these for each lipid individually, depending on the number of lipids and the number of beads of each lipid. You may need to change these variables (NMOL and NAT, respectively) in the script. You can also adapt the script to take data from only a subset of the lipids, for example to calculate error bars on your distributions. Now, start the script or take a look at the files first and then start the script.
$ [ gedit/vi ] anglelist
$ [ gedit/vi ] angledistr.sh
$ ./angledistr.sh
You will see that a new directory has been created: ANGLES contains all data of all the angles of all the lipids as a function of time (cga_n.m.xvg), numbering the angles (n) and lipids (m) from 0 onward. Angle 0 corresponds to the first angle in the file anglelist. The distribution of the pooled data for each angle in the file anglelist can be found in the file ANGLES/cga_n-distr.xvg, where n is a number. You may inspect these files using a plot program, e.g. xmgrace. If all is well, the four angle distributions should look like those given in the files cfMappedAngleDistributions.agr, and cfMappedAngleDistributions.png.
These are the target distributions that we want to obtain with the CG model.
3. CG simulation
As a starting point, go to the directory refine/take0. We will use the Martini parameters as defined in martini_v2.0_DIPC.itp, i.e., the angles for the angles around the double-bond beads are defined with the minimum energy at 120 degrees. Carry out a short CG simulation (starting from cg_mdstart.gro).
$ grompp -p cg.top -c cg_mdstart.gro -n run.ndx -f martini_md.mdp -o cg.tpr
$ mdrun -deffnm cg
Alternatively, you may want an educated guess based on the distributions before you start out on take0. There is a tool to help you. Go back to the directory refine. The python script do_fit.py fits harmonic or periodic potentials to target distributions and returns parameters for the reference length or angle and force constant. NOTE: you will need to have the python package symfit installed! As an example, try the first angle:
$ python do_fit.py -f ANGLES/cga_0-distr.xvg -x0 100 -k 25 --harmonic --radians
You will see a window with the target distribution in blue and the first guess gaussian distribution in red. On the bottom of the window are sliders; you can adjust them to give a better fit. The parameter values will serve as the starting values of a best fit. Closing this window triggers the fitting procedure to start, and will result in best fit parameters and a new window showing the fit. If you remember to specify the --radians flag on the command line, the force constant should have the correct Gromacs units.
4. Calculate CG distributions and 5. iterate
By now, this is now easy! All you need is the same input (anglelist) and script (angledistr.sh) as for the mapped trajectory, except to change the name of the CG trajectory (to cg.xtc). Copy them over, edit, and analyze. After comparing the angle distributions of interest, adjust the relevant parameters in the CG .itp file and run a new simulation. It is probably a good idea to do every simulation in a new directory. Comparing distributions can be done using xmgrace, for example in directory refine/take0 after analysis:
$ xmgrace ANGLES/cga_0-distr.xvg ../ANGLES/cga_0-distr.xvg
A python tool that gives you a graphical comparison and quantitative comparison is also available:
$ python compare.py ANGLES/cga_0-distr.xvg ../ANGLES/cga_0-distr.xvg
Tools and scripts used in this tutorial
gromacs (http://www.gromacs.org/)
insane.py (downloadable here)
backward.py (downloadable here)
MDAnalysis package (pypi.python.org/pypi/MDAnalysis/0.8.1)
[1] Marrink, S. J., De Vries, A. H., and Mark, A. E. (2004) Coarse grained model for semiquantitative lipid simulations. J. Phys. Chem. B 108, 750–760.
[2] Marrink, S. J., Risselada, H. J., Yefimov, S., Tieleman, D. P., and De Vries, A. H. (2007) The MARTINI force field: coarse grained model for biomolecular simulations. J. Phys. Chem. B 111, 7812–7824.
[3] Balgavy, P., Dubnicková, M., Kucerka, N., Kiselev, M. A., Yaradaikin, S. P., and Uhrikova, D. (2001) Bilayer thickness and lipid interface area in unilamellar extruded 1,2-diacylphosphatidylcholine liposomes: a small-angle neutron scattering study. Biochim. Biophys. Acta 1512, 40–52.
[4] Risselada, H. J., and Marrink, S. J. (2008) The molecular face of lipid rafts in model membranes. Proc. Natl. Acad. Sci. U.S.A. 105, 17367