This page will show the procedure to perform Replica Exchange Simulation with MMTSB toolset.
As I mentioned earlier, REX can be performed in vacuum and implicit solvent, if you want you can run in explicit solvent, but the speed will be very slow and not realistic to do so, and more importantly, mmtsb did not consider explicit exchange when they developed the toolset.
The first step is to prepare the starting structures, the starting structures can be either native structure or random extended structures. For native structure, download from PDB protein data bank, and formatted by command:
minCHARMM.pl -par cmap,blocked 1E0Q.pdb > 1e0q.pdb
remember to check the energy of the minimized structure, for native structure this is not a problem, but for generated extended structure, the default settings may not enough, you need to perform longer minimization by setting minsteps and sdsteps. Some time it could be unreasonable structure, and the energy would be very high, in this case, you need regenerate the structures.
To generate extended structures, run script:
genini.sh 1e0q.one 24 [explicit|anything else]
24 is # of copies, explicit will also generated solvated proteins, if not give any letter besides explicit, or the cshell will exist since there is no arguments.
If explicit were given, after the structures were generated, 1e0q.[0..9].wb.pdb could be deleted.
Check Grymoire to see the pitfall of Cshell and BASH, it will be confusing switching between different shells, so make a decision and stick to it.
To check the energy of the structures, run in tcsh:
set i 1
while ($ <= 24)
enerCHARMM.pl -par cmap,blocked 1e0q.$i.pdb
@ i++
end
After the starting structures were generated, we will be able to run REX simulation with command:
aarex.pl -n 1000 -par archive,natpdb=1e0q.pdb -temp 24:300:850 -log rex.log -charmmlog charmm.log -elog energy.log -mdpar cmap,blocked,(gb) -dir gbrex_va 1e0q.{1,2,3,.......24}.pdb
Once the simulation complete, the first thing is to create DCD trajectory files by:
trajarchive2dcd.pl [options] -dir directory options: [-ref refpdbfile] options: [-inx :]
This command will create dcd trajectories under each replica (./aa*), if -ref was given, it will create aa*-ref-rmsd.dat and aa*-ref-qscore.dat under each replica.
To load the trajectories into VMD, run Command:
vmd -e ~/toolset/loadrex.tcl -args gbrex_va [drawcell]
After the loading of the trajectories, open tcl script window, run following command:
display resetview
mol fix all
mol free 0
translate to -0.5 -0.5 0
scale by 0.35
then use mouse put the molecule into the corresponding cell, then repeat the commands.
(1)Making movie
Use the ppmtompeg, with the parameter file, in this file, the lower the value of IQ PQ BQ, the better the quality of movie, I usually set to 1, and if needed, use ffmpeg to convert to mp4, the command is
ffmpeg -i *.mpg *.mp4
(2) RMSD GRID Plot
There are many ways to plot the rmsd vs timesteps, rmsd grid is the most neat way, run script:
rmsd_grid.sh gbrex_va
This will create rmsd_grid.png in ./gbrex_va/.
(3) Energy Overlap Plot
If we want to look the energy distribution of all replicas or temperatures, run script:
ener_overlap.sh [-byclient|-bytemp] [dir of replica]
This will create ener_overlap_replica.png or ener_overlap_temp.png in ./gbrex_va/
(4) 2D WHAM Free Energy Analysis
Theory about the PMF and Gibbs Free Energy Calculation see the page of Free Energy Calculation.
MMTSB has the power to do WHAM analysis, in order to do WHAM analysis, reaction coordinates were needed for each temperature at each replica exchange cycle, let's define the first reaction coordinate is radius of gyratio and second is native contacts.
following are the contents of rgyr.function:
sub analyze {
my $mol=shift;
return &Analyze::radiusOfGyration($mol,1);
}
1;
following are the contents of qscore.function:
sub analyze {
my $mol=shift;
my $refmol=Molecule::new("1e0q.pdb"); #need to change the name of reference structure accordingly
my $analyze=Analyze::new($refmol);
my $qsc=$analyze->qscore($mol);
return $qsc->{all};
}
1;
then run command:
rexanalysis.pl -dir gbrex_va -inx 1:1000 -function rgyr.function > rgyr.data
Once we have the data for the reaction coordinates, we can run WHAM analysis by:
rexanalysis.pl -dir gbrex_va -inx 1:1000 -wham rgyr:rgyr.data:5:15:20=qscore:qscore.data:0:1:20 -whamtemp 300:400:500:600:700:800
PMF files will be created containing two-dimensional PMF profiles at each temperature.
the profiles can be plotted with gnuplot, here we show a method to plot the graph and make a movie.
run wham_multi.sh to create wham.gnu, and then run gnuplot wham.gnu, following are the contents in wham_multi.sh:
#!/bin/tcsh
set i=300
echo "reset" > wham.gnu
echo "set terminal gif animate" >>wham.gnu
echo "set output 'rgyr_qscore.gif'" >>wham.gnu
echo "set zrange [0:150]" >>wham.gnu
echo "set cbrange [0:150]" >> wham.gnu
while ($i <= 800)
echo "set dgrid3d 50,50,1" >>wham.gnu
echo "set pm3d">>wham.gnu
echo "set contour base">>wham.gnu
echo "splot 'rgyr_end2end.$i.pmf' u 1:2:3 w pm3d notitle">>wham.gnu
echo "unset dgrid3d" >>wham.gnu
echo "unset pm3d" >>wham.gnu
echo "unset contour">>wham.gnu
echo "\!sleep 1">>wham.gnu
@ i+=100
end
(5)Umbrella replica exchange run
To enhance the sampling of conformations, we add biased potential to the energy function, as was described in the theory of REX in parent page. this bias could be anything you are interested, mostly it is three dimensional constraints, let's first take a look at distance constraint.
1E0Q is a small hairpin, we use residue distance as constraint, the distance from 5A to 50A of Ca of first residue to Ca of last residue, the condition file looks like this:
bias resd
300 force=1.0,target=5,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=17,atom2=CA
300 force=1.0,target=6,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=17,atom2=CA
............
300 force=1.0,target=49,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=17,atom2=CA
300 force=1.0,target=50,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=17,atom2=CA
The first has to be assigned to what kind of bias function it is, followed by condition for each replica., they are all in 300K, and the force constant is 1.0 kcal/mol. Umbrella replica exchange use aarex.pl with same options except assign -condifile fileofcondition:
aarex.pl -n 1000 -condfile cond_resd -par archive,natpdb=1e0q.pdb -log rex.log -charmmlog charmm.log -elog energy.log -mdpar cmap,blocked,(gb) -dir gbrex_va_resd 1e0q.pdb
Once the simulation is complete, we can use WHAM tool to analyze the data. here I use Alan Grossfield's wham code. in order to do the analysis, we need to generate a number of input files:
rexinfor.pl -inx 1:1000 -condsel 0 run,val1 -dir gbrex_va_resd >data.r5
run this command for remaining conditions. In addition to these files we also need a master input file for WHAM analysis. It contains the names of data files and information of the midpoints and force constants of the umbrella potential for each condition:
data.r5 5.0 1.0
data.r6 6.0 1.0
..........
data.r49 49.0 1.0
data.r50 50.0 1.0
we can perform this by tcsh:
set i=0
set j=0
while ($j <= 50)
rexinfo.pl -inx 1:1000 -condsel $i run,val1 -dir gbrex_va_resd >data.r$j
@ i ++
@ j++
end
let's name this file whaminput_resd, then run command:
wham 5 50 60 0.00001 300 0 whaminput_resd whamoutput_resd
The options are 5A to 50A, 60 binis for the histogram, a convergence threshold 0.00001, temperature 300K, "0" meaningless but required,and input and output files.
Plot the output file to see the free energy as a function of end_to_end distance.
It is also possible to generate two dimensional PMFs from one dimensional umbrella run by introduction a second reaction coordinate and carrying out 2D WHAM analysis, but here we perform 2D umbrella replica exchange simulation.
here we use dihetral angle of 7.CA,8.CA,9.CA,10.CA as second reaction coordinate,
bias resd
bias dihe
300 force=1.0,target=5,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=17,atom2=CA force=0.05,target=-150,sel1=7.CA,sel2=8.CA,sel3=9.CA,sel4=10.CA
300 force=1.0,target=5,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=17,atom2=CA force=0.05,target=-100,sel1=7.CA,sel2=8.CA,sel3=9.CA,sel4=10.CA
300 force=1.0,target=5,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=17,atom2=CA force=0.05,target=-50,sel1=7.CA,sel2=8.CA,sel3=9.CA,sel4=10.CA
300 force=1.0,target=5,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=17,atom2=CA force=0.05,target=0,sel1=7.CA,sel2=8.CA,sel3=9.CA,sel4=10.CA
300 force=1.0,target=5,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=17,atom2=CA force=0.05,target=50,sel1=7.CA,sel2=8.CA,sel3=9.CA,sel4=10.CA
300 force=1.0,target=5,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=17,atom2=CA force=0.05,target=100150,sel1=7.CA,sel2=8.CA,sel3=9.CA,sel4=10.CA
300 force=1.0,target=5,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=17,atom2=CA force=0.05,target=150,sel1=7.CA,sel2=8.CA,sel3=9.CA,sel4=10.CA
.....................
300 force=1.0,target=50,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=17,atom2=CA force=0.05,target=150,sel1=7.CA,sel2=8.CA,sel3=9.CA,sel4=10.CA
Once the simulation has completed, we can extract both reaction coordinate from the rexserver.data with command:
rexinfo.pl inx 1:1000 condsel 0 run,val1,val2 dir rex2d > data.r5d-150
And we also need to create input file whaminput_2dx which has the following:
data.r5d-150 5.0 150.0 1.0 0.000015
data.r5d-100 5.0 100.0 1.0 0.000015
data.r5d-50 5.0 50.0 1.0 0.000015
data.r5d0 5.0 0.0 1.0 0.000015
data.r5d50 5.0 50.0 1.0 0.000015
data.r5d100 5.0 100.0 1.0 0.000015
data.r5d150 5.0 150.0 1.0 0.000015
............
data.r50d100 6.0 100.0 1.0 0.000015
data.r50d150 6.0 50.0 1.0 0.000015
The small second force constant for dihedral term is due to the convertion from radians to degrees 0.05 by (Pi*Pi)/(180*180).
Then run command:
wham-2d Px=0 5 50 60 Py=360 -180 180 60 0.00001 300 0 whaminput-2d whamoutput-2d
Plot the PMF, compare to 1D PMF.