Advanced Search
Last updated date: Nov 4, 2020 Views: 2458 Forks: 0
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−4bar−1.
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−1nm−2). 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
Related files
Category
Do you have any questions about this protocol?
Post your question to gather feedback from the community. We will also invite the authors of this article to respond.
Tips for asking effective questions
+ Description
Write a detailed description. Include all information that will help others answer your question including experimental processes, conditions, and relevant images.
Share
Bluesky
X
Copy link