non-CA atoms in nma

[Xinqiu Yao]

I received an Email asking a question about NMA. I post it below and happy to hear your comments. Thanks!

---

Hello,

I was wandering if there was a way to analyse non-CA atoms in the NM analysis; in particular, I am trying to include the porphyrin ring of a pigment. If so, I'll make sure to get the response onto the bio3d forum. I am just thinking it might not be possible.

Thanks!

PS

script:

library(bio3d)
# Define function for spring force constants
"my.ff" <- function(r, ...) {
  ifelse( r>10, 0,r^(-2) )
}

pdb <- read.pdb("fmo1.pdb")

# Indices for CA atoms
a.inds<-atom.select(pdb,elety = "CA")#, chain = c("A","D","E"), operator = ("AND"))
b.inds<-atom.select(pdb,elety = c("C1B","C2B","C3B","C4B","CHB","C1A","C2A","C3A","C4A","CHA","CBD","CAD","C3D","C4D","C2E","C1D","CHD","C1C","C2C","C3C","C4C","CHC"))
z.inds<-a.inds
z.inds$atom <- c(a.inds$atom,b.inds$atom)
z.inds$xyz <- c(a.inds$xyz,b.inds$xyz)

# Build hessian matrix
h <- build.hessian(pdb$xyz[ ,z.inds$xyz ], aa.mass , pfc.fun=my.ff)

# Diagonalize and obtain eigenvectors and eigenvalues
modes <- eigen(h, symmetric=TRUE)

[Lars Skjærven]

This looks good. The tricky business will probably be the parametrization of reasonable spring force constants for such a model. i.e. is it reasonable to assume that you can use the same force field for the ligand atoms as for the CA atoms? Note that build.hessian doesn't do the mass weighting anymore, so you can omit the aa.mass argument if you have a new version of bio3d. This was moved to a separate function in 2.1 I think. Also, you might find useful a set of internal functions in nma.pdb.R for producing nma objects.

We have also a experimental version of all-atom NMA. See branch feature_aanma and function aanma().

[Alexander Fokas]

Hello Lars,

For starters I think it would be fine to use the same force field. However, I do know that force fields do exist for the pigments (https://simtk.org/home/psii_ff), so I guess that would be something to consider after.

Will I be unable to investigate the frequencies if there are no mass weightings? This was my ultimate goal. If not, should the above script get me there?

The all atom looks very interesting, I'll be sure to have a look once this is sorted out.

Thanks, Alexander

[Xinqiu Yao]

Hi Alex,

Yes, you are right. Without mass you won't get frequencies but force constants (So called "energetic modes"). I think it is still better to use nma.pdb() but modify your pdb file by converting the non-CA atoms that you would like to include to "CA" atoms. For example,

spdb <- trim(pdb, inds=z.inds)
write.pdb(spdb, file="enm.pdb")

# open the file enm.pdb and change all non-CA atoms to "CA" manually (resid can be arbitrary, e.g. "LIG")

work.pdb <- read.pdb("enm.pdb")

# replace "XX" to mass value you want. For example 12 if you are representing every carbon atom in the ring.
modes <- nma.pdb(work.pdb, pfc.fun=my.ff, mass.custom=list("LIG"=XX))

Lars, do you have any comment for this? Thanks!

[Lars Skjærven]

Perhaps.. you can also do it step wise as suggested with internal functions.

library(devtools)
load_all("~/workspace/bio3d/ver_devel/bio3d")

# read protein, and strip all but CA and ligand
pdb <- read.pdb("1rx2")
sele <- atom.select(pdb, "calpha", resid="FOL", operator="OR")
pdb <- trim(pdb, sele)

# define force field
ff <- load.enmff("anm")

# build hessian
hess <- build.hessian(pdb$xyz, ff)

# mass weight hessian
masses <- atom2mass(pdb)
hess_mw <- .nma.mwhessian(hess, masses=masses)

# diagonalize mass weighted hessian
ei <- .nma.diag(hess_mw)

# produce final nma object
m <- .nma.finalize(ei, xyz=pdb$xyz, temp=300, masses=masses,
                   natoms=nrow(pdb$atom), keep=200, call=NULL)

Note the need of sourcing the nma.pdb.R file here (done through the devtools package in this case). You can also just copy the functions from here. I'm not sure why these functions are not exported.. Xinqiu?

[Lars Skjærven]

or, if you want, fetch the effective hessian of the CA atoms:

# NMA on protein with ligand. outputs on CA modes.
pdb <- read.pdb("1rx2")
sele <- atom.select(pdb, "calpha", resid="FOL", operator="OR")
pdb.in <- trim(pdb, sele)

# define force field
ff <- load.enmff("anm")

# output modes for CA atoms only
out.inds <- atom.select(pdb.in, "calpha")
pdb.out <- trim(pdb.in, out.inds)

# fetch number of atoms and sequence
natoms.in  <- ncol(pdb.in$xyz)/3
natoms.out <- ncol(pdb.out$xyz)/3
sequ <- pdb.in$atom$resid

# fetch masses
masses.in <- atom2mass(pdb.in)
masses.out <- atom2mass(pdb.out)

# build hessian, extract effective hessian
hessian <- .nma.hess(pdb.in$xyz, init=list(pfcfun=ff),
                     hessian=NULL, inc.inds=out.inds)

# mass weight
hessian <- .nma.mwhessian(hessian, masses=masses.out)

# diagonalize
ei <- .nma.diag(hessian)

# make the NMA object
m <- .nma.finalize(ei, xyz=pdb.out$xyz, temp=300, masses=masses.out,
                        natoms=natoms.out, keep=200, call=NULL)


# NMA without the ligand present
m2 <- nma(pdb, ff="anm")

# compare modes
rmsip(m, m2)

does it make sense?

[Xinqiu Yao]

Hi Lars,

All internal functions starting with '.' (period) will not be exported. This is set on the top of "NAMESPACE" file. I am not sure if this is designed so at beginning, but I guess it would be good to have "internal functions" internally used. Right?

[Barry Grant]

That was our intention - we can add exceptions if necessary though...

[Alexander Fokas]

Great, thanks for your help!