Extracting High-Frequency Motions From the NMFF Pathway
Integrating the low-resolution models with atomically detailed simulations and subsequently determining thermodynamic quantities (e.g., free energy) along the pathway is of extreme interest in the field of computational biophysics. Here, we present a way to integrate information from low-resolution models with the atomically detailed models.
Specifically, we will perform minimization with restraints starting from each structure along the pathway generated using NMFF. We will use CHARMM for this purpose. The example CHARMM input file (kinase_bead1.inp) that acomplish's this task is as follows:
* This tutorial example illustrates how to perform minimization with restraints
* in CHARMM.
* Note first we will build the missing atoms in the starting structure.
* In an accompanying example we will do all this with mmtsb command.
! Read in the parameter and topology files.
! We'll use the $CHARMMDATADIR as our source for top/param files
open unit 1 read form name "$CHARMMDATA/top_all22_prot_cmap.inp"
read rtf card unit 1
close unit 1
open unit 1 read form name "$CHARMMDATA/par_all22_prot_cmap.inp"
read param card unit 1
close unit 1
! Read the sequence explicitly
!read sequ card unit 5
!* Sequence for 214 residue ADK
read sequ card
* Adenylate Kinase
MET ARG ILE ILE LEU LEU GLY ALA PRO GLY ALA GLY LYS
GLY THR GLN ALA GLN PHE ILE MET GLU LYS TYR GLY ILE
PRO GLN ILE SER THR GLY ASP MET LEU ARG ALA ALA VAL
LYS SER GLY SER GLU LEU GLY LYS GLN ALA LYS ASP ILE
MET ASP ALA GLY LYS LEU VAL THR ASP GLU LEU VAL ILE
ALA LEU VAL LYS GLU ARG ILE ALA GLN GLU ASP CYS ARG
ASN GLY PHE LEU LEU ASP GLY PHE PRO ARG THR ILE PRO
GLN ALA ASP ALA MET LYS GLU ALA GLY ILE ASN VAL ASP
TYR VAL LEU GLU PHE ASP VAL PRO ASP GLU LEU ILE VAL
ASP ARG ILE VAL GLY ARG ARG VAL HSD ALA PRO SER GLY
ARG VAL TYR HSD VAL LYS PHE ASN PRO PRO LYS VAL GLU
GLY LYS ASP ASP VAL THR GLY GLU GLU LEU THR THR ARG
LYS ASP ASP GLN GLU GLU THR VAL ARG LYS ARG LEU VAL
GLU TYR HSD GLN MET THR ALA PRO LEU ILE GLY TYR TYR
SER LYS GLU ALA GLU ALA GLY ASN THR LYS TYR ALA LYS
VAL ASP GLY THR LYS PRO VAL ALA GLU VAL ARG ALA ASP
LEU GLU LYS ILE LEU GLY
generate PRO0 first NTER last CTER setup
open unit 9 read card name adk_mov1.pdb ! Read coordinates of the first structure
read coor pdb unit 9
close unit 9
! Since the structures from the NMFF path do not have hydrogens atoms
! we add them here. Use ic build and then hbuild.
coor init select hydrogen end
! For consistency with defaults in the mmtsb tool minCHARMM.pl, we use
! this hbuild set-up
hbuild atom cdie eps 80.0 cutnb 10.0 ctofnb 7.5 ctonnb 6.5 shift vshift bygr
update atom rdie eps 1 cutnb 12 ctofnb 11 ctonnb 8 switch vswitch bygr
! Now we will do a short minimization with harmonic restraints on CA atoms
!Restraint CA atoms with high force constant
! Copy current coordinates to comparison set
coor copy compare
cons harm force 100 mass select ( ( resid 1:214 .and. segid PRO0 ) ) .and. ( atom * * ca ) end
mini sd nprint 10 step 0.02 nstep 500 ! Steepest-Descent minimization
mini abnr nprint 10 step 0.005 nstep 1000 tolenr 1e-05 ! ABNR minimization scheme
open unit 10 write form name "movmin1.pdb" ! minimized coordinates
write coor pdb unit 10
close unit 10
Use the MMTSB tool convpdb.pl to "convert" from generic resiude naming convention to "charmm22" convention for residues and atoms.
convpdb.pl -out charm22 mov1.pdb > adk_mov1.pdb
The ADK conformation is minimized using CHARMM with the command:
$CHARMMEXEC < kinase_bead1.inp > kinase_bead1.out
Preferably, use minCHARMM.pl routine in the MMTSB Tool Set to perform the same task.
minCHARMM.pl -par minsteps=0,sdsteps=10,sdstepsz=0.02 \
-par minsteps=20,minmode=abnr \
-par trunc=switch,dielec=rdie,cutnb=12,cuton=8,cutoff=11 \
-par xtop=top_all22_prot_cmap.inp \
-par xpar=par_all22_prot_cmap.inp \
-cons ca self 1:214_100 -cmd mmtsbmin.inp -log mmtsbmin$i.log \
mov$i.pdb > movmin$i.pdb
Now minimize other intermediate states along the pathway using the same protocol. You have a choice of either submitting jobs manually after making changes to CHARMM input file (kinase_bead1.inp) to read the next structure along the pathway or execute a simple shell script (runcharmm.sh) to automate this task. If multiple CPU's are available all structures can be minimized in parallel. This is particularly attractive feature of the method and makes it computationally tractable to investigate conformational transitions of large biomolecules. Note that here we performed very short minimiation to save time but in practince one should run longer minimization till gradient of about 1e-06 or lower is achieved.
After the minimization jobs are finished, load in VMD all the minimized intermediate conformations along the ADK pathway. You can do this by issuing the following commands in the "tkcon" console of VMD:
animatepdbs 1 30 "movmin%0d.pdb"
Turn on the simultaneously the representation of the molecule to "Cartoon" for the entire protein and to "Licorice" for only few selected residues, numbered 36, 88, 123, 156, 16 and 13 in the PDB files. Now play the trajectory and visualize the atomically detailed residue motions that were missing in your initial NMFF pathway.
In principle, it is also feasible to extract thermodynamically relevant information from the NMFF pathways. Although, to put this into practice much more computationally intensive calculations are required and is beyond the scope of this workshop. If time permits, interested participants may work on the suggested exercise.