Molecular Dynamics simulations

Molecular Dynamics simulations

For readers new to computer simulations, we suggest the reader introductory GROMACS

http://www.mdtutorials.com/gmx/

and MARTINI

http://www.cgmartini.nl/index.php/tutorials-general-introduction-gmx5

tutorials.

In this protocol, we assume that the reader has already installed the desired version of GROMACS simulation package. In our article, we have performed all molecular dynamics simulations using GROMACS 5.0.5

Peptide preparation

The MARTINI model cannot describe protein folding due to the missing possibility to form backbone hydrogen bonds. Thus, we have to define the protein secondary structure a priori and this structure is kept throughout the whole simulation. These were the performed steps:

1) Obtain an atomistic structure

We can obtain the initial peptide structure by several means, for example:

a) Download from the PDB database.

wget http://files.rcsb.org/download/2MLT.pdb

NOTE: The file can contain additional molecules that should be removed before proceeding further.

b) Using the MODELLER program.

https://salilab.org/modeller/

For alpha-helical peptide, following script can be used:

https://salilab.org/modeller/wiki/Make%20alpha%20helix

c) Using the Builder GUI in PyMOL molecular viewer.

2) Convert the atomistic structure to coarse-grained representation

Subsequently, use python "martinize" script to generate the peptide topology and structure files based on the atomistic model. The script is provided by the authors of the MARTINI force field:

http://cgmartini.nl/index.php/tools2/proteins-and-bilayers/204-martinize

If the experimental structure is available, use the DSSP program to determine the peptide secondary structure automatically:

python martinize.py -f peptide.pdb -o topol.top -x peptide_CG.pdb -dssp /full/path/to/dssp -p backbone -ff martini22

Alternatively, specify custom secondary structure:

python martinize.py -f peptide.pdb -o topol.top -x peptide_CG.pdb -ss "file_or_string" -p backbone -ff martini22

Use "HHHHHHHHHCHHHHHHHHH" string to assign alpha-helical structure with fully flexible kink in the center to a 19-residue long peptide.

NOTE: Unnatural residues or alternative residue names (e.g., HID, HIE) may not be recognized correctly.

3) Energy minimization in vacuum

Place the peptide into a cubic box (separated by at least 1 nm from the boundaries):

            gmx editconf -f peptide_CG.pdb -d 1.0 -o peptide_CG.gro

Perform the energy minimization using the steepest descent algorithm (integrator = steep) with convergence criterion set to 10.0 kJ mol−1nm−1 (emtol = 10):

gmx grompp -p topol.top -f min.mdp -c peptide_CG.gro -o peptide_min.tpr

gmx mdrun -v -deffnm peptide_min

Membrane preparation

We can conveniently obtain assembled lipid bilayer from CHARMM-GUI website (both topology and structure files for GROMACS are provided). By convention, we assemble the lipid bilayer in the XY-plane. Since peptides commonly have non-zero net charge, we can generate the system without ions or remove the solvent altogether.

System preparation

1) System assembly - configuration

Place peptides equally on both membrane leaflets (in-plane orientation at the level of lipid headgroups) at the desired peptide-to-lipid ratio. Rotate the peptides to orient their hydrophobic regions towards the hydrophobic core. For the individual placement of peptides, use:

gmx editconf -f peptide_CG.pdb -o peptide_N.gro -box <box> -center <pos> -c -rotate <rot> -princ

, where <box> are the membrane-system dimension, <pos> is the position of peptide center of mass, and <rot> is the peptide rotation.

Finally, merge all configuration files into system_vac.gro and edit the second line of the gro file to state the correct number of particles in the system.

Note: First line of the gro file is a comment, second line is the number of particles in the system, and the last line are the box dimensions.

2) System assembly - topology

Now, we have a structure file and we have to make changes to the topology to match it. Make the following changes to the membrane topology file:

a) Include the peptide-parameter itp file in the file header.

b) Write the number of peptides in the [ molecules ] (reflecting also the position within the structure file).

c) Energy minimization in vacuum

Then we remove the steric clashes between the added peptides and lipids:

gmx grompp -p topol.top -f min.mdp -c system_vac.gro -o system_vac_min.tpr

gmx mdrun -v -deffnm system_vac_min

Note: Check for possibility of threaded molecules through ring groups. The energy minimization will finish, but the subsequent simulation will crash.

3) System assembly - solvation

Next step is the addition of solvent into the energy-minimized peptide-membrane system. To get the correct system density, we will need to adjust the default particle size (-radius parameter). If necessary, we can prevent the insertion of water molecules into the hydrophobic core by increasing the size of selected particles.

Copy the vdwradii.dat file (found in the GROMACS installation path) to the working directory and adjust the vdW radii. Increase the vdW radii (up to 1 nm) for both membrane and peptide particles to fill all voids inside the membrane.

From the MARTINI force field website, download a box of pre-equilibrated MARTINI water (water.gro) and run:

gmx solvate -cp system_vac_min.gro -cs water.gro -o system_sol.gro -radius 0.21 -p topol.top

Note: The topology topol.top file should now also list the correct number of water molecules.

4) System assembly - adding ions

To match our experimental conditions, we add 100 mM of NaCl. We also add excess of ions to neutralize system net charge.

Firstly, create an empty file empty.mdp since no simulation parameters are needed for the addition of ions. Then, create a tpr file that is necessary to add ions to the solvated system:

gmx grompp -f empty.mdp -c system_sol.gro -p topol.top -o add_ions.tpr

Finally, add the ions by replacing randomly selected water molecules with ions:

gmx genion -s add_ions.tpr -o system_sol_ion.gro -p topol.top -pname NA -nname CL -conc 0.100 -neutral

5) Energy minimization in solvent

Perform energy minimization just like in previous steps:

gmx grompp -p topol.top -f min.mdp -c system_sol_ion.gro -o system.tpr

gmx mdrun -v -deffnm system

System equilibration

Use leap-frog algorithm for integrating Newton’s equations of motion. Generate the initial particle velocities (gen_vel = yes) for the temperature of 323 K (gen_temp = 323) according to Maxwell distribution. To keep the temperature constant, use velocity rescaling thermostat (modified with a stochastic term). For proper temperature distribution, use two coupling groups for peptide-membrane and solvent particles. To generate the index groups, run:

gmx make_ndx -f system.gro -o index.ndx << EOF

Protein | POPC

W | Ion

EOF

Maintain system pressure (an membrane tension) via Berendsen thermostat (ref_p = 1.0).

For systems with lipid bilayers, use semi-isotropic coupling scheme. Use Verlet cutoff scheme set the vdW radius to 1.1 nm. Set compressibility of the simulation box to 3.10⁻⁴bar⁻¹.

Perform the equilibration in five subsequent steps with different simulation times:

(1) 0.5 ns  (time step of 2 fs)

(2) 1.25 ns (time step of 5 fs)

(3) 1 ns    (time step of 10 fs)

(4) 30 ns   (time step of 20 fs)

(5) 15 ns   (time step of 20 fs)

In steps (1)--(4), apply postional restraints on the peptide backbone beads (force constant of 1000 kJ mol⁻¹nm⁻²). Due to the changes of the box size, turn on the scaling of the restraint reference positions (refcoord_scaling = com). In step (5), perform a short production molecular dynamics was run without any positions restraints.

For the first step, run:

gmx grompp -f equi-1.mdp -n index.ndx -c system.gro -p topol.top -o npt-1.tpr -maxwarn 1

gmx mdrun -rdd 1.6 -v -deffnm npt-${i}

For subsequent runs, include the checkpoint file in the grompp command using the parameter -t.

Production run

Use the final configuration (.gro) checkpoint (.cpt) to start the production simulation run.

In the config (.mdp) file, change the barostat to Parrinello-Rahman with coupling time constant 12 ps. Then run:

gmx grompp -f md.mdp -n index.ndx -c npt-5.gro -p topol.top -o md.tpr

gmx mdrun -rdd 1.6 -v -deffnm md

Category