1.2 All-atom MD simulation of POPC lipid bilayers

Files for this tutorial (charmm-gui.gtz)
Optional input file to wrap atoms into unit cell (tutorial1.2-crd_convert.inp)

Cell membranes have many biological roles such as substrate transport, signal transduction, and cell adhesion. By using GENESIS, the user can simulate biological membranes at the atomistic level. One of the difficulties in such simulations is initial setup of the bilayer structure or equilibration of the system. Because each lipid molecule has a random structure, it may take much time to equilibrate the whole system, if the initial structure has an ordered conformation. CHARMM-GUI is one of the useful tools to generate initial structures of lipid bilayers for MD simulations (S. Jo et al., 2007 ). It provides the automatic membrane builder for not only pure lipid bilayers but also protein-membrane systems. Recently, GENESIS format has been introduced into the CHARMM-GUI membrane builder as well as NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM formats (J. Lee et al., 2015 ). Here, we explain how to simulate POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) lipid bilayers with GENESIS and CHARMM-GUI using the CHARMM C36 force fields (J. Klauda et al., 2010 ). Other sample control files for membrane simulations are also available here.

1.2.1 Setup initial structure with CHARMM-GUI

We use the Bilayer builder in the CHARMM-GUI. If the server is not working, files prepared by us included in the tutorial files (available from the link just below the title of this page) can be used instead. In this case, this section can be skipped. Now, let us build a POPC lipid bilayer. At the first step of the bilayer builder, we choose “Membrane Only System”:


In the next step, we decide the number of lipids and water molecules to be added in the system. Here, we put 40 POPC and 40 POPC lipids in the upper and lower leaflets, respectively (Note that these numbers might be a lower limit for stable MD simulation with spdyn). Because each lipid molecule has a certain area per lipid (Aexp = 68.3Å2 for POPC), the initial box size (X, Y) = (52.27Å, 52.27Å) can be estimated from sqrt{(number of lipids in one leaflet)*(area per lipid)}. tutorial1.2-step2

In the next steps, we add 0.15 M KCl in the system, mimicking the cellular environment, with the Monte Carlo method for ion placing. Then, we click “Next Step” three times. We check “GENESIS” in “Input Generation Options” and select the NPT ensemble at T = 303.15K in “Equilibration Options”. At this temperature, POPC bilayers show liquid phase rather than gel phase (Tm = 271K). After clicking the “Next Step”, we can download a set of input files for minimization, equilibration, and production runs for GENESIS at the final step of the membrane builder.


1.2.2 Simulation with GENESIS

First of all, let us see the PDB structure in the downloaded file by using VMD:

# make a new directory to do this tutorial and copy the downloaded file
$ mkdir tutorial-1.2
$ cd tutorial-1.2
$ cp ~/Download/charmm-gui.tgz ./

# check the downloaded files
$ tar -xzf charmm-gui.tgz
$ cd charmm-gui
$ vmd step5_assembly.pdb


We use this PDB file as the initial coordinates of the simulation. In the current working directory, there is the genesis directory, where several control files for GENESIS are contained:

$ cd genesis
$ ls
restraints                 step6.2_equilibration.inp  step6.6_equilibration.inp
step5_charmm2genesis.str   step6.3_equilibration.inp  step7_production.inp
step6.0_minimization.inp   step6.4_equilibration.inp
step6.1_equilibration.inp  step6.5_equilibration.inp

These control files are mainly composed of three steps: energy minimization (step6.0), equilibration (step6.1-6.6), and production runs (step7). In total, 8 simulations are carried out step by step. Basically, we do not need to edit these control files (If the user uses GENESIS1.0, pme_ngrid_x, y, and z should be set in the [ENERGY] section). During Step6.0-6.6, dihedral angle restraints are applied on some parts of the system to avoid unexpected structural change in the lipid molecule such as isomerization or cis-trans transition of double bonds. In the restraints directory, such restraint functions are defined, and they are called from the control files. Positional restraints are also imposed on z-coordinates of lipid phosphorous atoms in order to prevent undulation of the membrane plane. These restraint forces are gradually reduced to zero during the equilibration.

Let us check parameters in the control files (step6.0-step7) before running the simulation. For step7, we can see:

$ less step7_production.inp
forcefield        = CHARMM    # [CHARMM]
electrostatic     = PME       # [CUTOFF,PME]
switchdist        = 10.0      # switch distance
cutoffdist        = 12.0      # cutoff distance
pairlistdist      = 13.5      # pair-list distance
pme_nspline       = 4
water_model       = NONE
vdw_force_switch  = YES

integrator        = VVER      # [LEAP,VVER]
timestep          = 0.002     # timestep (ps)
nsteps            = 500000    # number of MD steps
crdout_period     = 1000
eneout_period     = 1000      # energy output period
rstout_period     = 1000 
nbupdate_period   = 10
rigid_bond        = YES       # constraints all bonds involving hydrogen
fast_water        = YES 
shake_tolerance   = 1.0D-10
ensemble          = NPT       # [NVE,NVT,NPT]
tpcontrol         = LANGEVIN  # thermostat and barostat
temperature       = 303.15
pressure          = 1.0       # atm
gamma_t           = 1.0
gamma_p           = 0.05
isotropy          = SEMI-ISO
type              = PBC       # [PBC]

All control files already include optimal or well-established MD parameters. For example, in the [ENERGY] section, we specify switchdist = 10 Å, cutoffdist = 12 Å, and vdw_force_switch = YES, which are recommended in the membrane simulations with the CHARMM C36 force fields (J. Lee et al., 2015 ). We use the SHAKE/SETTLE algorithms for rigid bonds, and Langevin thermostat and barostat for temperature and pressure control. As for the pressure control, semi-isotropic pressure coupling is used. The equations of motion are integrated by the velocity Verlet method.

Let us carry out energy minimization, equilibration, and production runs with spdyn. Here, we use 16 MPI processors as an example. You can use a smaller or larger number of MPI processors as well. Since the system size is relatively large, it may take more than 10 minutes in some equilibration runs and more than 3 hours in the production run.

# 1000 step energy minimization
# with dihedral angle restraint (k = 250 kcal/mol)
#      z-positional restraint on P atoms (k = 2.5 kcal/mol/Å2)
$ export OMP_NUM_THREADS=1
$ mpirun -np 16 /home/user/genesis/bin/spdyn step6.0_minimization.inp > log6.0

# 25ps MD in the NVT ensemble
# with dihedral angle restraint (k = 250 kcal/mol)
#      z-positional restraint on P atoms (k = 2.5 kcal/mol/Å2)
$ mpirun -np 16 /home/user/genesis/bin/spdyn step6.1_equilibration.inp > log6.1

# 25ps MD in the NVT ensemble
# with dihedral angle restraint (k = 100 kcal/mol)
#      z-positional restraint on P atoms (k = 2.5 kcal/mol/Å2)
$ mpirun -np 16 /home/user/genesis/bin/spdyn step6.2_equilibration.inp > log6.2

# 25ps MD in the NVT ensemble
# with dihedral angle restraint (k = 50 kcal/mol)
#      z-positional restraint on P atoms (k = 1.0 kcal/mol/Å2)
$ mpirun -np 16 /home/user/genesis/bin/spdyn step6.3_equilibration.inp > log6.3

# 100ps MD in the NPT ensemble
# with dihedral angle restraint (k = 50 kcal/mol)
#      z-positional restraint on P atoms (k = 0.5 kcal/mol/Å2)
$ mpirun -np 16 /home/user/genesis/bin/spdyn step6.4_equilibration.inp > log6.4

# 100ps MD in the NPT ensemble
# with dihedral angle restraint (k = 25 kcal/mol)
#      z-positional restraint on P atoms (k = 0.1 kcal/mol/Å2)
$ mpirun -np 16 /home/user/genesis/bin/spdyn step6.5_equilibration.inp > log6.5

# 100ps MD in the NPT ensemble without any restraints
$ mpirun -np 16 /home/user/genesis/bin/spdyn step6.6_equilibration.inp > log6.6

# 1000ps MD in the NPT ensemble
$ mpirun -np 16 /home/user/genesis/bin/spdyn step7_production.inp > log7

At the last step, we obtain the coordinates trajectory file step7_procution.dcd.

1.2.3 Visualization of the trajectory

Let us visualize the trajectory file by using VMD:

$ vmd ../step5_assembly.pdb -psf ../step5_assembly.psf -dcd step7_production.dcd

In the simulation, we can see that hydrophobic lipid tail groups fluctuate significantly and form random conformation. Water molecules are interacting with hydrophilic lipid head groups through hydrogen bonds. One of the important features of the lipid bilayers is that they are almost impermeable to ions because of the low dielectric environment in the hydrophobic core region.

In the VMD viewer window, we may find H-H bond in the TIP3P water, but it does not matter because it comes from definition of H-H bond in the PSF file for the SHAKE/SETTLE algorithms. We can also see that water molecules are spread out from the simulation box, but they are actually existing in the periodic image cells. To wrap molecules into the unit cell, we utilize crd_convert in the GENESIS analysis tool set. First, we prepare a control file for crd_convert (tutorial1.2-crd_convert.inp  ; same as the file in the top of this page), where pbc_correct = MOLECULE is specified in the [OPTION] section:

check_only     = NO          # (YES/NO)
trjout_format  = DCD         # (PDB/DCD)
trjout_type    = COOR+BOX    # (COOR/COOR+BOX)
trjout_atom    = 1           # atom group
pbc_correct    = MOLECULE    # (NO/MOLECULE)

and then execute crd_convert by the following command:

$ /home/user/genesis/bin/crd_convert tutorial1.2-crd_convert.inp

We obtain a converted trajectory file (step7_production_conv.dcd).

1.2.4 Trajectory analysis

We analyze structural properties of lipid bilayers. Among them, the area per lipid is one of the important quantities. For pure lipid bilayer systems, it can be calculated by Box_size_x*Box_size_y/(number of lipid molecules). Here, we have 40 POPC lipids in one leaflet, and Box_size_x and Box_size_y are output in the column 17 and 18 in the “INFO:” line of the log file:

$ less log7
            RMSG            BOND           ANGLE    UREY-BRADLEY        DIHEDRAL
        IMPROPER         VDWAALS           ELECT     TEMPERATURE          VOLUME
            BOXX            BOXY            BOXZ          VIRIAL        PRESSURE
         PRESSXX         PRESSYY         PRESSZZ
 --------------- --------------- --------------- --------------- ---------------
INFO:       1000          2.0000     -17944.5309     -31100.6338      13156.1030
         13.9359       1118.7362       4598.0045       1814.1727       4417.2590
         53.2784       -512.8975     -42589.1871        302.1313     190320.8275
         50.3044         50.3044         75.2098      -8696.4821         26.7517
        145.9381       -362.8987        297.2158

In order to plot the time courses of the area per lipid, we first run the grep command to search the string “INFO:” in the log file and create a new file containing those lines, and then execute gnuplot:

$ grep "INFO:" log7 > energy.log
$ gnuplot
gnuplot> set xlabel 'Time (ps)'
gnuplot> set ylabel 'Area per lipid (A^2)'
gnuplot> plot 'energy.log' using 3:($17*$18/40) with lines


 If your plot is different from our results, it may be due to different random number used in the simulations.

Although the calculated area per lipid is slightly smaller than the experimental value (Aexp = 68.3Å2 for POPC), GENESIS reproduces the averaged area per lipid of 64.7Å2 reported in the CHARMM C36 paper (J. Klauda et al., 2010 ). Longer MD simulation (~100 ns) is needed to obtain reliable structural data of lipid bilayers.


Written by Takaharu Mori@RIKEN Theoretical molecular science laboratory
June 30, 2016