Common applications of the MMTSB toolset

Page currently under construction!!!

Preparing protein simulations (from PDBfile to CHARMM trajectory):

This tutorial illustrates the step-by-step instructions to perform explicit solvent molecular dynamics of a protein solvated in a cubic box of water at 0.15 mM ionic NaCl concentration using the MMTSB toolset. The set of files and scripts provided below can be easily adapted to perform dynamics of any other protein of interest with different solvent box shapes and sizes. Here, we use RNAase peptide as an example.

Step-by-step instructions are as follows:

1. Download RNAse A crystal structure, 1RNU.pdb from the Protein DataBank.

2. Extract first 13 residues using convpdb.pl, centered and with segnames and in charmm22 format:

            convpdb.pl -center -sel 1:13 -segnames -out charmm22 -nohetero 1RNU.pdb > 1rnu_cpep.pdb

3. Visualize peptide fragment with VMD to make sure its what you want.

4. Minimize

            minCHARMM.pl -par minsteps=0,sdsteps=100,sdstepsz=0.02 \
          -par trunc=switch,dielec=rdie,cutnb=12,cuton=8,cutoff=11 \
          -par param=22x,cmap\
          -par xtop=top_all27_prot_na.rtf \
          -par xpar=par_all27_prot_na.prm \
          -par blocked,nter=none,cter=ct3 \
          -cons heavy self 1:13_5 -log mmtsbmin.log \
            1rnu_cpep.pdb > 1rnu_cpep_mmtsbmin.pdb

5. Write protein structure file (PSF) for use in later steps involving: (1) viewing structures in VMD and (2) determining how many ions to add for making the charge on the system neutral.

                genPSF.pl  -par param=22x,cmap\
              -par xtop=top_all27_prot_na.rtf \
              -par xpar=par_all27_prot_na.prm \
              -par blocked,nter=none,cter=ct3 \
               1rnu_cpep_mmtsbmin.pdb > 1rnu_cpep_mmtsbmin.psf

In the following steps 6-8 of the tutorial, we will use the MMTSB tools to solvate our minimized peptide. First we'll add solvent and ions to the system. Then we'll minimize the complete system with harmonic restraints on the peptide. The system will then be ready for molecular dynamics runs.

6. Add solvent box

First use convpdb.pl to solvate and add segment names - extract boxsize!

                convpdb.pl -out charmm22 -solvate -cubic -cutoff 8 1rnu_cpep_mmtsbmin.pdb | \
          convpdb.pl -out charmm22 -segnames > 1rnu_cpep_solv.pdb       

In this case, the boxsize is 38.743553 Å cube. Note down this number.

7. Add ions at 0.15 M concentration to the cubic box of water generated above

Counterions are added to solvated system by specifying the number of positive (SOD) and/or negative ions (CLA). Determine how many and which type of ions we need to neutralize this system using the CHARMM script how_many_ions_to_add.inp.

            $CHARMMEXEC pdbfile=1rnu_cpep_mmtsbmin psffile=1rnu_cpep_mmtsbmin \
          boxsize=38.74355 <how_many_ions_to_add.inp> how_many_ions_to_add.log &

Note: Use PSF file of the protein only and NOT one with the solvent at this step.

Now grep information about the number of ions from the log file "how_many_ions_to_add.log"

Number of positive ions (SOD) to add:

            grep "NPOS ->"  how_many_ions_to_add.log | tail -1

Number of negative ions (CLA) to add:

            grep "NNEG ->"  how_many_ions_to_add.log | tail -1

Now finally add ions to the system with waterbox (in this case 5 SOD and 6 CLA ions):

            convpdb.pl -ions SOD:5=CLA:6 1rnu_cpep_solv.pdb > 1rnu_cpep_solvions.pdb

Generate complete PSF of the system including waterbox and ions.

             genPSF.pl  -par param=22x,cmap,\
            -par xtop=top_all27_prot_na.rtf \
            -par xpar=par_all27_prot_na.prm \
            -par blocked,nter=none,cter=ct3 \
              1rnu_cpep_solvions.pdb>1rnu_cpep_solvions.psf

8. Now minimize with restraints and SHAKE constraints using PME electrostatics

             minCHARMM.pl -par minsteps=0,sdsteps=500,sdstepsz=0.02 \ 
             -par trunc=switch,cutnb=12,cuton=8,cutoff=11 \ 
             -par param=22x,cmap\ 
             -par xtop=top_all27_prot_na.rtf \ 
             -par xpar=par_all27_prot_na.prm \ 
             -par blocked,nter=none,cter=ct3 \ 
             -par shake,boxsize=38.743553,nblisttype=bycb \ 
             -cons heavy self 1:13_5 \ 
             -cmd mmtsbminsolvate.inp -log mmtsbminsolvate.log \ 
              1rnu_cpep_solvions.pdb > 1rnu_cpep_solvions_min.pdb

9. Run dynamics on this peptide in water using the mmtsb tool mdCHARMM.pl.

In this section we will take the configuration we just generated and run 2 ps of molecular dynamics on this system:

             mdCHARMM.pl  -par dynsteps=1000,dynens=NPT,dynitemp=298,dyneqfrq=1000 \
            -par dynnose=1,dynoutfrq=10,dynpress=1,echeck=20000 \ 
            -par trunc=switch,cutnb=12,cuton=8,cutoff=11 \
            -par param=22x,cmap \
            -par xtop=top_all27_prot_na.rtf \ 
            -par xpar=par_all27_prot_na.prm \
            -par blocked,nter=none,cter=ct3 \
            -par shake,boxsize=38.743553,nblisttype=bycb \
            -cmd mmtsbdynsolvate.inp -log mmtsbdynsolvate.log \
            -enerout 1rnu_cpep_d0.ene -trajout 1rnu_cpep_d0.dcd \
            -restout 1rnu_cpep_d0.res -final 1rnu_cpep_d0.pdb \
             1rnu_cpep_solvions_min.pdb

10. Restart dynamics from previous step and perform production dynamics using NVT ensemble:

             mdCHARMM.pl  -par dynsteps=1000,dyntemp=298,\ 
            -restart 1rnu_cpep_d0.res -final - 1rnu_cpep_solvions_min.pdb \
            -par trunc=switch,cutnb=12,cuton=8,cutoff=11 \
            -par param=22x,cmap\
            -par xtop=top_all27_prot_na.rtf \
            -par xpar=par_all27_prot_na.prm \
            -par blocked,nter=none,cter=ct3 \
            -par shake,boxsize=38.743553,nblisttype=bycb \
           -cmd mmtsbdynsolvate1.inp -log mmtsbdynsolvate1.log \
            -enerout 1rnu_cpep_d1.ene -trajout 1rnu_cpep_d1.dcd \
            -restout 1rnu_cpep_d1.res

11. Visualize the dynamics trajectory.

After running the dynamics, visualize the trajectory with vmd. To look at molecular dynamics trajectories with vmd, you need to read in a pdb and a dcd file. The following command will accomplish this.

            vmd 1rnu_cpep_solvions_min.pdb 1rnu_cpep_d1.dcd

12. Analyze the trajectory

Process dcd trajectory to extract root mean square deviation of the peptide CA atoms with respect to the reference PDB after least-squares superposition as a function of time.

            analyzeCHARMM.pl -rms -sel CA -comp 1rnu_cpep_solvmin.pdb -psf \
            1rnu_cpep_solvmin.psf 1rnu_cpep_d1.dcd >rms.out & 

Visualize data using Gnuplot software.

            gnuplot
            pl ‘rms.out’ u 2:3 w lp

13. For further analysis of the trajectory, see the MMTSB Tool Set - Trajectory Analysis tutorial.

Preparing protein simulations for replica-exchange:

replace text

Visualizing the electrostatic surface potential of a macromolecule:

First, we select the protein from a PDB file (1enh.pdb), add hydrogen atoms, and center the molecule at the origin. Then, we calculate the electrostatic potential on a grid and generate the molecular surface of the protein. Finally, we view the molecule in VMD, projecting the electrostatic potential onto the molecular surface.

            convpdb.pl -nsel protein 1ENH.pdb  |  complete.pl  |  convpdb.pl -center > 1enh.center.pdb
            pbCHARMM.pl -emap phi.dx 1enh.center.pdb
            pbCHARMM.pl -dx -grid epsx grid.dx 1enh.center.pdb
            vmd 1enh.center.pdb phi.dx grid.dx

In vmd, select:

            Graphics/Representations/Drawing Method [Surf]
            Graphics/Representations/Coloring Method [Volume]

For a more thorough description, see the MMTSB Tool Set - Continuum electrostatics calculations tutorial.

For more examples:

Download and follow tutorials that were prepared for past MMTSB workshops.