iRED analysis is described in Prompers, J. J., Bruchweiler, R.

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

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

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

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`

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.

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)

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`

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