Protocol Sheet: ACETK.BLDAMBER.v1

Build a globular protein in the Amber forcefield

Software requirements: VMD.1.8.7, AmberTools 1.5
Knowledge requirements: Use of VMD and Amber LEaP
Input requirements: We provide all input documents in the attachment, including the PDB of the structure we are going to use. As long as the software mentioned above is installed properly, you should be ready to go. Please note that this is not a tutorial on how to parameterize a ligand, new residue, etc. Please see ACETK.PRMAMBER for instructions on how to do that.

Download ACETK.BLDAMBER.v1

Disclaimer. This protocol is based on a best practice. In no way, we guarantee that it is optimal for your system. You use it at your own risk.

Suggested directory structure:

acetk.bldcharmm.v1
├── build
│   ├── 1LKK.pdb
│   ├── build-functions.tcl
│   ├── determine_ions.vmd
│   ├── par_files
│   │   ├── Y2P.frcmod
│   │   └── Y2P.off
│   ├── sleap_1lkk.in
│   ├── sleap_1lkk_ionize.in
│   └── split.vmd
├── equil
└── run

Suggestion. After building use protocol ACETK.EQ to equilibrate your system.

STEP1: Prepare the PDB file(s)

When beginning to build any system, you may start with one or multiple PDB files, depending on what you are trying to do. While in this tutorial we are only dealing with one starting PDB file, this procedure can apply to multiple input structures.

1) Split the system(s) into desired components. We want to remove the specific components that we want and discard anything superfluous. For the protein and ligand, we will strip the hydrogens and let the Amber software guess the coordinates based on the parameter files. We will also preserve the crystal waters from the PDB structure, and this is something that should be done unless there is a clear reason not to (even if the positions of the hydrogens are unknown).

 vmd -dispdev text -e split.vmd

2) Perform additional modifications, where necessary. In our case, because the naming conventions in the PDB and in our parameter file for the ligand are different, we must do some manual work to make sure they match up. Amber LEaP will give an error otherwise.

sed -E 's/PTR/Y2P/' ligand.pdb | \
sed -E 's/\ OH\ /\ OG\ /'      | \
grep -v 'C   ACE B 256'        | \
awk '$6 ~ /256/ {gsub(/256/, "251")}; 1' > ligand_fixed.pdb   

# Command 1: Change name of Phosphotyrosine residue to match frcmod file.
# Command 2: Change name of an oxygen, again to match that of frcmod file.
# Command 3: Remove troublesome carbon atom and allow leap to geuss coordinates.
# Command 4: Change residue 256 to 251 since it is part of that residue.

Important: We are dealing with a very peculiar pdb file with the ligand. The atoms of the ACE residue are split into two "residues", and the atoms are not listed sequentially. This can happen when the crystalographer sees an atom in two positions, or just can't determine the exact position. We have removed atom 46 with the script above after visual inspection. Furthermore, sleap requires the atoms of a residue to be grouped together (even with non sequential atom numbers), otherwise it is incapable of dealing with this issue. That means that you have to manually move atom 47 from the bottom of the file and put it after the second atom. sleap will throw an error otherwise. Always be aware that the software you are working with may make any number of assumptions that are difficult to imagine, so check your work carefully for errors.

STEP2: Preliminary Build

Run sleap in order to build a preliminary system. This will include the salvation box, ion neutralization, and sulfide bond formation. Amber already has all standard residues parameterized in it's force field. However, additional residues or ligands must be parameterized on their own (see ACETK.PRMAMBER) and the parameter files must be loaded into sleap in the script before the PDB is read. We have provided the complete script, but we highly suggest looking it over to see what each command does. It can then easily serve as a template for future systems.

Run sleap:

sleap -f sleap_1lkk.in

STEP3: Determine Ion concentration

One peculiarity of Amber LEaP is that it does not provide a way to automatically ionize a system to a desired salt concentration or ionic strength. It can automatically neutralize the system, but as yet it can not add ions to a desired salt concentration. Therefore, if you want to do simulations at physiologically accurate salt concentrations, you must manually calculate and add those ions. I have made a script that manually calculates monovalent salt concentration based on the water volume of the system. It also counts how many salt ions were added during the first run of sleap so that you can figure out how many you will have to adjust for.

To run the script, type the following command.

vmd -dispdev text -e determine_ions.vmd

STEP4: Ionization and Final Structure

Now we want to add the salt ions. For this, we have to run sleap all over again. We have supplied a different input file for this, so that you do not have to modify the other one. Use your favorite text editor and edit sleap_1lkk_ionize.in so that 31 sodium ions are added and 27 chlorine ions are added. The addions part of your input file should look like this:

...
addions fullmol Na+ 31
addions fullmol Cl- 27
...

Why 31 and 27? We saw from running our script in step three that we needed 29 sodium ions and 29 chlorine ions solely based on the water concentration. But we also already had 4 sodium ions added in the first build by sleap so as to neutralize the system. So we must subtract those evenly from both ions so that we can maintain the charge of the 4+ added to balance out negative charges on the protein, all while keeping an accurate ionic strength. So we subtract 2 from both the recommended sodium and chlorine totals, and add the four to the number of sodiums to add, giving 31 and 27. What would happen if we had an odd number of sodium atoms balancing our system, such as 3? Having balanced charges is the first concern, so we would simply choose a number of ions to add that would give 3+ ions while either just under- or over-shooting the correct ionic strenght (e.g. 30+/27- or 31+/28-). Simulations are finite systems, so we simply have to make the best compromise possible.

Now we are ready for the final step. Run the command below, and the system is finally fully built.

sleap -f sleap_1lkk_ionize.in

General note: sleap is extremely picky about what imperfections it can accept as input, and is still in development by the AmberMD team. In the past it has done such strange things as ask for the angle between two atoms (Atom1-Atom2-Atom1). tleap, the predecessor of sleap, is much less so. Therefore, if you have problems with sleap, one potential solution is simply to use tleap. The input files are usually interchangeable. However, being older software, it is possible that tleap is missing errors that sleap has been designed to catch. Always check your finished system manually and watch for any unexpected behavior during simulations.