Experiment of Rplica Exchange Simulation

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++


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);

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};

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: 

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



(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++


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.r50d­100 6.0 ­100.0 1.0 0.000015
data.r50d­150 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.