AMHGoModelAndUmbrellaSampling

Nicholas Schafer edited this page Nov 17, 2017 · 5 revisions

Introduction

This page introduces the use of AMH-Go, a non-additve structure based model included in the AWSEM simulation package, and the use of umbrella sampling over a collective coordinate for the purposes of performing free energy profile and expectation value calculations. AMH-Go is a non-additive structure based model, which is described in detail in: Eastwood, Michael P., and Peter G. Wolynes. "Role of explicitly cooperative interactions in protein folding funnels: a simulation study." The Journal of Chemical Physics 114 (2001): 4702.

Example: running an AMH-Go simulation

Obtain a structure of your protein of interest from the PDB (www.pdb.org) or elsewhere.

The structure should have coordinates for at least the backbone atoms and C-beta for all residues that you are interested in simulating. For the purposes of illustration, we will use Ubiquitin (pdb id: 1UBQ). Download the pdb file to a local directory. Use some tools, or edit the pdb file using a text editor, so that it only includes the chain that you want to simulate. In this case we will be simulating chain A, the only chain in this particular pdb file. If you are following along in this tutorial, rename this modified file to 1ubqA.pdb.

Generate the input files and copy the necessary parameter files to the current working directory.

Next, run the AWSEM script to start generating the input files:

PdbCoords2Lammps.sh 1ubqA 1ubqA

The first argument is the name of the pdb file that you generated and modified but without the ".pdb" extension. The second argument will serve as the pre/postfix for the input files that are generated by the script. In this case, we have chosen to use the same name for both, but this is not necessary. Note that many of the scripts, if run with no arguments, will show the basic syntax, such as:

$ PdbCoords2Lammps.sh

> /home/ns24/opt/script/PdbCoords2Lammps.sh pdb_id project_name

The files generated include coord, data, seq and in files. The coord file is a temporary file used only for constructing the data file and may be safely deleted. The in file is the LAMMPS input script. The seq file contains the one letter amino acid sequence, and the data file contains generic information about the number of atoms and bonds in the simulation, the number of atom and bond types, the dimensions of the simulation box, the masses of the atoms, the atom IDS, chain IDs, residue IDs, atom types, atom charges and x, y and z coordinates of the atoms as well as bond coefficients and a list of bonds between atoms.

Next, copy amh-go.gamma and fix_backbone_coeff.data from the parameters directory to your local directory. Edit the fix_backbone_coeff.data file so that it contains only the following blocks, and make sure that none of the blocks have a "-" after the closing bracket, as this will inactivate the corresponding potential:

[Chain]
10.0 10.0 30.0
2.45798 2.50665 2.44973

[Chi]
20.0 -0.83

[Epsilon]
1.0

[Rama]
2.0
5
 1.3149  15.398 0.15   1.74 0.65 -2.138
1.32016 49.0521 0.25  1.265 0.45  0.318
 1.0264 49.0954 0.65 -1.041 0.25  -0.78
    2.0   419.0  1.0  0.995  1.0  0.820
    2.0  15.398  1.0   2.25  1.0  -2.16

[Rama_P]
3
 0.0    0.0 1.0   0.0  1.0   0.0
2.17 105.52 1.0 1.153 0.15  -2.4
2.15 109.09 1.0  0.95 0.15 0.218
 0.0    0.0 1.0   0.0  2.0   0.0
 0.0    0.0 1.0   0.0  2.0   0.0

[SSWeight]
0 0 0 1 1 0
0 0 0 0 0 0

[ABC]
0.483 0.703 -0.186
0.444 0.235 0.321
0.841 0.893 -0.734

[AMH-Go]
1.0
2.5
8.0

In the [AMH-Go] block, the first parameter is the scaling of this term relative to the other terms. The second parameter is the "p-value", or non-additivity exponent. Higher p-values give high barriers and more cooperative folding. Values in the range of 2-3 have been found to give realistic levels of cooperativity. The third and final parameter is the cutoff for defining native contacts between residues, and is in units of Angstroms. For the purposes of this tutorial, we will use the values p=2.5 and a cutoff of 8 angstroms.

The amh-go.gamma file contains information about the strength of interactions between residues of different identities and sequence separation classes. The default file is consistent with those parameters used in the original study where the model was introduced and subsequent studies. The model is a "vanilla" model is the sense that the strength of interaction between residues does not depend on their identity. This can be easily changed by modifying the amhgo-gamma file, but that is beyond the immediate scope of this tutorial.

We now need to generate the amh-go.gro file, which will be used to determine which residues are in contact in the native state. Make sure that you renumber the PDB file that you will be using to generate the amh-go.gro file such that the first residue index is 1. Then, the file can be generated using:

python /path/to/awsemmd/frag_mem_tools/Pdb2Gro.py 1ubqA amh-go.gro

where, again, this assumes that you have a file called 1ubqA.pdb and creates a file callled amh-go.gro. The gro file contains the number of atoms, residue IDs, residue types, atom types, atom IDs, and x, y and z coordinates of the atoms.

Finally, it is necesary to generate an "ssweight" file, which contains information about how to bias the secondary structure using the Ramachandran potentials. Since we don't want to provide any extra bias to the native structure in this case, we can just use a file containing 76 (the number of residues in 1UBQ chain A) of "0.0 0.0", meaning that the extra bias to an alpha or beta conformation are both zero for all residues.

$ for i in {1..76}; do echo "0.0 0.0"; done > ssweight

Run a simulation.

You are now ready to run your first simulation with the AMH-Go model. By issuing the command

$ lmp_serial < 1ubqA.in

you should see output that looks something like

LAMMPS (9 Oct 2012)
Scanning data file ...
  3 = max bonds/atom
Reading data file ...
  orthogonal box = (-10000 -10000 -10000) to (10000 10000 10000)
  1 by 1 by 1 MPI processor grid
  228 atoms
  302 bonds
Finding 1-2 1-3 1-4 neighbors ...
  5 = max # of 1-2 neighbors
  5 = max # of special neighbors
76 atoms in group alpha_carbons
76 atoms in group beta_atoms
76 atoms in group oxygens
Chain flag on
Chi flag on
Rama flag on
Rama_P flag on
SSWeight flag on
ABC flag on
AMH-Go flag on

Setting up run ...
Memory usage per processor = 0.891418 Mbytes
Step Temp E_pair E_mol TotEng Press Volume 
       0          300    69.685911    2804.7907    2656.0812    -19470.21    37318.056 
    1000     309.6291    3.2541839    154.84377    -57.57167    377.89729    37318.056 
    2000    317.92045    4.5022927    145.15843   -46.009096    -513.3088    34061.896 
    3000    289.04461    4.2000368    136.29229    -80.36671   -53.054637    34061.896 
    4000    284.42523    3.7390817    122.59595   -81.326705   -701.07713    32570.479 
    5000    338.64666    4.4512755    122.41484   -38.536759     533.1893     35143.81 
    6000    284.88008    1.6497933    105.21846   -101.70821    708.35035    37609.352 
    7000    309.85035    2.2742725    107.32425   -71.231702      319.249    36253.327 
    8000    288.93468    4.5076069    114.54767   -69.546086     144.3742    37828.582 
    9000    318.63316    3.7868783    101.90289   -60.801079   -147.73095    37141.499 
   10000    317.62052    2.9372909    106.25488   -54.292209    172.91441    36277.457 
Loop time of 5.3855 on 1 procs for 10000 steps with 228 atoms

Pair  time (%) = 0.740053 (13.7416)
Bond  time (%) = 0.0497167 (0.923159)
Neigh time (%) = 0.00521493 (0.0968328)
Comm  time (%) = 0.00129962 (0.0241318)
Outpt time (%) = 0.00178814 (0.0332029)
Other time (%) = 4.58743 (85.1811)

Nlocal:    228 ave 228 max 228 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost:    0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs:    10891 ave 10891 max 10891 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs:  21782 ave 21782 max 21782 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 21782
Ave neighs/atom = 95.5351
Ave special neighs/atom = 2.64912
Neighbor list builds = 10
Dangerous builds = 0

Example: applying a collective coordinate bias

Add a fix qbias command to your input file.

In order to sample different parts of the landscape, it is often useful to apply biases to your simulation Hamiltonian. In this case, we will apply a bias to the Q coordinate. Low values of Q correspond to unfolded structures, while high values of Q correspond to folded, near-native structures. To apply a bias to your molecule, you can add code like this to your input file:

# apply a qbias
fix	      qbias alpha_carbons qbias fix_qbias_coeff.data # apply a Q-biasing potential as specified in fix_qbias_coeff.data
fix_modify    qbias energy no # do not add the energy from the Q-biasing potential to the total potential energy
variable      biasinge equal f_qbias

See the LAMMPS manual for general information about how "fixes" work. In this case, the first line applies a fix with the label "qbias" to the "alpha_carbons" group, which is of type "qbias" and takes as input the file "fix_qbias_coeff.data". The qbias fix assumes that you have a file named "rnative.dat" in the current directory which contains a square matrix of native CA-CA distances. Qbias fix uses that matrix to calculate the instantaneous Q value of the simulated structure, the biasing energy, and the biasing forces at every timestep. The procedure for generating rnative.dat is described below. The second line, starting with "fix_modify", prevents the potential energy coming from the bias from being automatically added to the total potential energy for the purposes of output. Whether or not you want to automatically include the biasing energy into the total potential energy depends on how you will be subtracting out the bias later. A complete explanation of this is beyond the scope of this tutorial, but most WHAM-like algorithms explicitly calculate the biasing energy internally and add it to the unbiased total potential energy, so including it here would be tantamount to double-counting. The last line creates a variable so that the biasing energy can be output separately if so desired.

The fix_qbias_coeff.data file has a format similar to that of fix_backbone_coeff.data discussed previously. An example fix_qbias_coeff.data file looks like:

[QBias_Exp]
KQBIAS
Q0
2
0.15

In this case KQBIAS should be replaced by a floating point number that dictates the strength of the biasing potential. For models with larger barriers it is necessary to apply larger biasing potentials, but applying too larger of a biasing potentials will reduce the overlap between adjacent sampling "umbrellas" and thereby impact the quality of the unbiasing calculation. In general it is necessary to check whether or not your sampling is sufficient by (at least) checking for significant overlap between sampling umbrellas. Values of KQBIAS between 100 and 1000 are appropriate for different models/proteins/purposes when using AWSEM. The Q0 parameter should replaced by a number between 0 and 1 and specifies the center of the biasing potential. The next two lines indicate that the biasing potential is harmonic (an exponent of 2) and that the width of the gaussian well used for biasing increases with increasing sequence separation according to |i-j|^0.15.

Compute collective coordinates on the fly using compute functions

AWSEM and LAMMPS can be used to compute some collective coordinates on-the-fly during the simulation using LAMMPS' compute functionality. An example block of input code that computes several collective coordinates including a Q value, radius of gyration, contact Q, the number of total contacts, and the total potential energy, is given below (replace SAMPLETIME with how often you want to output to the wham.dat file, e.g., 1000):

# output wham.dat, a file with collective variables that can be used to compute pmfs
compute     qw alpha_carbons qwolynes rnative.dat 2 0.15
variable    qw equal c_qw
compute     rg alpha_carbons gyration
variable    rg equal c_rg
compute     qo alpha_carbons qonuchic cutoff 12.0 nativecoords.dat 1.2
variable    qo equal c_qo
compute     tc beta_atoms totalcontacts 6.5 2
variable    tc equal c_tc
variable    step equal step
variable    E_P equal pe

fix         wham all print SAMPLETIME "${step}        ${qw} ${rg} ${qo} ${tc} ${E_P}" file wham.dat screen no title "# timestep qw           rg          qo            tc            energy"

In order to use these compute functions, it is necessary to prepare the rnative.dat and nativecoords.dat files. This can be done using:

python /path/to/awsemmd/create_project_tools/GetCACADistancesFile.py 1ubqA rnative.dat
python /path/to/awsemmd/create_project_tools/GetCACoordinatesFromPDB.py 1ubqA nativecoords.dat

Note:Use GetCACADistanceFile_multiChain.py for multi-chain simulation. The result of running a simulation with code like that above in the input file should be a file called wham.dat in the directory where the simulation was run.

Setup multiple directories, Run all simulations, Collect the data

Once these extra items have been added to the input file, you are ready to run umbrella sampling by specifying an array of Q0 values over a set of directories. All of these simulations can be run in parallel, and once the simulations have finished you will want to collect the wham.dat files, and perhaps other files such as energy.log, in order to compute 1D and 2D free energy profiles and expectation value plots. Performing these calculations is beyond the scope of this tutorial, but you may find the following tools to be useful: WHAM.jar pyMBAR Grossfield WHAM

You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.
Press h to open a hovercard with more details.