Checking Loops RMSDs: MD vs OTF.
Rolando P. Hong Enriquez
25 June, 2015
1.- Overview
In this report we compared the loop that is crossed by the ligand (DIH) with respect to the other 2 loops in the trimeric protein PNP. Additionally a comparison with the same process in a MD simulation is also performed. All the rmsd values in this report are calculated using the same frame of reference, which is the first frame from the MD trajectory.
2.- Data
The data to be analyzed comes from two different sources:
A long MD trajectory that displays the binding of the ligand (DIH) to the binding pocket of the protein (PNP) using a path that we called gate01 in which the ligand crosses a protein loop before arribing to the binding site.
Images comming from the OTF string method calculation that represent converged binding paths of DIH to PNP. The number of images are 21 (from primary OTF calculations) and 40 (from OTF calculations derived from reparametrizations of the primary calculation). In these images, we also specify several force constants that were applied to the protein using ColVar module of namd to keep the protein from rotating/translating.
# get into the folder
wd0<-getwd()
wd1<-paste(wd0,'/RMSD_LOOPS', sep="")
setwd(wd1)
# load data
rmsdMD <- read.table("rmsdLoops_SergioMD.txt")
rmsd40K200 <- read.table("rmsdLoops_I40K200.txt")
rmsd40K500 <- read.table("rmsdLoops_I40K500.txt")
rmsd40K1000 <- read.table("rmsdLoops_I40K1000.txt")
colnames(rmsdMD) <- c("LoopA","LoopB","LoopC")
colnames(rmsd40K200) <- c("LoopA","LoopB","LoopC")
colnames(rmsd40K500) <- c("LoopA","LoopB","LoopC")
colnames(rmsd40K1000) <- c("LoopA","LoopB","LoopC")
head(rmsdMD,3)## LoopA LoopB LoopC
## 1 1.476330 1.411721 1.331597
## 2 0.000494 0.000503 0.000499
## 3 1.153333 1.456340 1.279710In these data files the column LoopA always represent the loop in PNP that is crosses by the ligand during the binding process.
3.- Visualization
#plot
#layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
par(mfrow=c(2,2))
boxplot(rmsdMD,ylab="RMSD",col=(c("red","gray","gray")),main="MD")
boxplot(rmsd40K200,ylab="RMSD",col=(c("red","gray","gray")),main="40K200")
boxplot(rmsd40K500,ylab="RMSD",col=(c("red","gray","gray")),main="40K500")
boxplot(rmsd40K1000,ylab="RMSD",col=(c("red","gray","gray")),main="40K1000")
A simple exploration of the data suggest that in all cases LoopA shows smaller RMSD from the reference structure, the other two loops have higher RMSD values and there is even a preserved order in RMSDs in all cases:
- RMSD(LoopA) < RMSD(LoopC) < RMSD(LoopB)
It is also to note that in general we also find higher RMSD values (for instance for LoopA) in the OTF calculations compared with the MD simulation.
#create data frame for LoopA in all conditions
maxLen<-length(rmsdMD$LoopA)
c1<-c(rmsd40K200$LoopA, rep(NA, maxLen - length(rmsd40K200$LoopA)))
c2<-c(rmsd40K500$LoopA, rep(NA, maxLen - length(rmsd40K500$LoopA)))
c3<-c(rmsd40K1000$LoopA, rep(NA, maxLen - length(rmsd40K1000$LoopA)))
all<-data.frame(rmsdMD$LoopA,c1,c2,c3)
colnames(all) <- c("MD","I40K200","I40K500","I40K1000")
boxplot(all,ylab="RMSD",col=(c("grey","blue","green","red")),main="LoopA")