iRED analysis

Instructions on how to perform iRED (isotropic reorientational eigenmode dynamic) analysis of bond vectors to obtain order parameters (S²). For Amber 10 or later.
iRED analysis is described in Prompers, J. J., Bruchweiler, R. J. Am. Chem. Soc. 2002,124, 4522-4534. DOI: 10.1021/ja012750u

The analysis is composed of two executable python scripts and a module that are located in /temp4/bio/Ired. The module can be extended such that other types of bond vectors can be handled. At the moment, backbone N-H, side-chain methyl, aromatic side-chains, and Bruchweiler-dictionary bond vectors are supported

The programs are in /temp4/bio/Ired

Samuel Genheden, 2011-2012

Skripten är fortfarande aktuella för att generera cpptraj input. Den sista raden i input-scripten bör ändras från
analyze matrix matired vecs 128 out bv.out
till 
analyze matrix matired order orderparamfile order.out

Det finns t ex kommandon för att räkna ut S2 direkt i cpptraj i stället för att post-processa ptraj output som mina skript gör. 
Sam 30/3-15

Prerequisite

It is assumed that you have at least one trajectory of the system of interest. If you want to perform several types of analysis on the protein (e.g. calculate different kinds of order parameters) it is recommended to strip the water molecules from the trajectories.

You need to create a reference pdb-file that is used to identify the bond vectors in your protein. It should have at least the same protein atoms and in the same order as the system in the trajectories. Therefore, it is recommended that it is created from the prmtop and prmcrd that was used to initiate the MD simulations. For instance:

changeparm << EOF prmtop p ref.pdb m prmcrd q EOF

Calculating the covariance matrix, eigenvalues and eigenvectors

Once the prerequisites are met, the covariance matrix of the bond vectors should be calculated. This, and the diagonalization of it to obtain the eigenvalues and eigenvectors are performed by ptraj. Amber10 or later should be used.

First,we need to create input to ptraj that tells which bond vectors that are of interest. This is created by a Python script called create_bv_inpt.py. Type,

create_bv_inpt.py nh ref.pdb bv.out mdcrd5 > ptraj.in

nh - the type of bond vectors of interest. See below, for the supported types of vectors.
ref.pdb - a reference pdb-file used to find the bond vectors. Created above
bv.out - a ptraj output-file containing the eigenvectors and the eigenvalues
mdcrd5 - the MD trajectory. You can specify one of more, separated by space. If you would like to only process a part of the trajectory, enclose the trajectory and the selection in quotation-marks, e.g., "mdcrd5 1 100"

ptraj.in now contains the following lines (or something similar):

trajin mdcrd5 vector v1 :2@N ired :2@H vector v2 :3@N ired :3@H ... matrix ired name matired order 2 analyze matrix matired vecs 128 out bv.out

The first line is a standard line that reads in the MD trajectory. Then follows 128 lines with bond vector definitions, where 128 is the total number of bond vectors. In the example above,we define backbone N-H bond vectors. The last two lines perform the actual creation of the covariance matrix and the calculation of the eigenvalues/eigenvectors, respectively. Here it is important that vecs are set to the total number of bond vectors.

Now, run ptraj

ptraj prmtop ptraj.in >& ptraj.out

Calculating the order parameters

Once ptraj is finished, the the eigenvalues and eigenvectors should be combined to calculate the order parameters. This is performed by a Python script called mat2s2.py. Type,

mat2s2.py nh ref.pdb offset bv.out

nh - the type of bond vectors of interest. See above

ref.pdb - a reference pdb-file. See above.

offset - a residue number offset, in case the real residue numbers do not start at 1. This number will be added to all residue numbers. Usually it is set to 0.

bv.out - a ptraj output-file containing the eigenvectors and the eigenvalues. See above. You can specify one or more of these if you would like to for instance average over windows.

The order parameters will be printed to standard output. If the more than one ptraj output-file was specified, the standard error will be printed as well.

Supported vector types

At the moment, the following vector types can be specified in create_bv_inpt.py and mat2s2.py:

nh - selects backbone N-H vectors
me - selects side-chain methyl vectors of Val, Thr, Leu, Ile, Met, and Ala residues
ar - selects N-H and C-H vectors of Phe, Tyr, Trp and His residues
ars - as ar, but the program will ask the user to select a single residue number
dic - selects side-chain vectors of the dictionary of Bruschweiler, described in Li, D-W., Bruschweiler, R. J. Am. Chem. Soc. 2009,131, 7226-7227. DOI: 10.1021/ja902477s
ala - selects the side-chain methyl of Ala (it is missing from the dictionary)

Extending the analysis scripts

Additional types of bond vectors can be introduced straightforwardly into the create_bv_inpt.py and mat2s2.py scripts by editing the def_bv.py module. This is a two-step process:

Create a function that is called find_bv_XXX.
XXX is a name of your choice. In fact, the function can be called almost anything, but it is recommended to call it like this.
This function should as its only input argument take a filename of a pdb-file.
This function should as its only output, return a list with all the bond vectors of interest. Each item in the list should contain in order, a) the residue name, b) the residue index, c) the first atom, and d) the second atom of the bond vector. All these variables should be strings.
Add the new functionality to the list of recognized bond-vector codes and to the list of bond-vector functions.
These defintions are at the end of the def_bv.py module. At the moment it looks like this:
bv_codes = ['nh','me','ar','ars','dic','ala'] bv_funcs = {} bv_funcs["nh"]=find_bv_nh bv_funcs["me"]=find_bv_me bv_funcs["ar"]=find_bv_ar bv_funcs["ars"]=find_bv_ars bv_funcs["dic"]=find_bv_dic bv_funcs["ala"]=find_bv_ala

bv_codes is a list of strings with the recognized bond-vector codes. Add XXX to this list.
bv_funcs is a map where the indices are the strings in bv_codes and the values are the functions.
Just add another line at the end, e.g.
bv_funcs["XXX"]=find_bv_XXX

Dictionary averaging

To calculate order parameters that are compatible with the Bruschweiler dictionary, the individual order parameters for a residue needs to be averaged. There exists a script to do this as well as calculate the dictionary entropies from the order parameters.

Execute

av_dic.py s2.out

where s2.out is the output from mat2s2.py on a set of order parameters calculated with dic as the bond-vector type. The script will write out the average order parameters as well as the entropy.

Sample sander.in files

To be used for any long simulation of proteins with PME and PBC
Use TIP4P water

sander.in1
Restrained minimization
&cntrl
irest=0,ntx=1,
imin=1,maxcyc=1000,drms=0.0001,ntmin=2,
ntc=2,ntf=2,
cut=8.0,
ntpr=100,ntwx=0,ntwv=0,ntwe=0,
ipol=0,igb=0,
ntr=1,restraint_wt=100,
restraintmask="!:WA= & !@H="
&end

sander.in2
Restrained NVT equilibration
&cntrl
irest=0,ntx=1,ig=-1,
nstlim=10000,dt=0.002,
temp0=300.0,ntt=3,gamma_ln=2.0,
ntc=2,ntf=2,
cut=8.0,
ntpr=500,ntwx=0,ntwv=0,ntwe=0,
ntb=1,
ipol=0,igb=0,
ntr=1,restraint_wt=50,
restraintmask="!:WA= & !@H="
&end

sander.in3
Restrained NPT equilibration
&cntrl
irest=1,ntx=5,
nstlim=10000,dt=0.002,
temp0=300.0,ntt=3,gamma_ln=2.0,
ntc=2,ntf=2,
cut=8.0,
ntpr=500,ntwx=0,ntwv=0,ntwe=0,
ntb=2,ntp=1,pres0=1.0,taup=1.0,
ipol=0,igb=0,
ntr=1,restraint_wt=50,
restraintmask="!:WA= & !@H="
&end

sander.in4
Equilibration, 1 ns NPT
&cntrl
irest=1,ntx=5,
nstlim=500000,dt=0.002,
temp0=300.0,ntt=3,gamma_ln=2.0,
ntc=2,ntf=2,
cut=8.0,
ntpr=500,ntwx=0,ntwv=0,ntwe=0,iwrap=1,
ntb=2,ntp=1,pres0=1.0,taup=1.0,
ipol=0,igb=0,
ntr=0
&end

sander.in5
Production 10 ns , sampling frequency 10 ps                                                            &cntrl
irest=1,ntx=5,
nstlim=5000000,dt=0.002,
temp0=300.0,ntt=3,gamma_ln=2.0,
ntc=2,ntf=2,
cut=8.0,
ntpr=5000,ntwx=5000,ntwv=0,ntwe=0,iwrap=1,
ntb=2,ntp=1,pres0=1.0,taup=1.0,
ipol=0,igb=0,
ntr=0
&end

sander.in4 if many independent simulations
Equilibration NPT with seed
&cntrl
irest=0,ntx=1,ig=28529,
nstlim=100000,dt=0.002,
temp0=300.0,ntt=3,gamma_ln=2.0,
ntc=2,ntf=2,
cut=8.0,
ntpr=500,ntwx=500,ntwv=0,ntwe=0,
ntb=2,ntp=1,pres0=1.0,taup=1.0,
ipol=0,igb=0,
ntr=0
&end