Simulating a pharmaceutical compound using antechamber and the Generalized Amber Force Field

By Ross Walker and Sishi Tang
Updated for AMBER 18 by Michael Barton and Tyler Luchko

Learning Outcomes

Introduction

Antechamber is designed to be used with the "general AMBER force field (GAFF)"¹. This force field has been specifically designed to cover most pharmaceutical molecules and is compatible with the traditional AMBER force fields in such a way that the two can be mixed during a simulation. Like the traditional AMBER force fields, GAFF uses a simple harmonic function form for bonds and angles but unlike the traditional protein and DNA orientated AMBER force fields the atom types used in GAFF are much more general such that they cover most of the organic chemical space. The current implementation of the GAFF force field consists of basic atom types and special atom types. The charge methods used can be HF/6-31G* RESP² or AM1-BCC³.

By design, GAFF, is a complete force field (so that missing parameters rarely occur), it covers almost all the organic chemical space that is made up of C, N, O, S, P, H, F, Cl, Br and I. Moreover, since GAFF is totally compatible with the AMBER macromolecular force fields it should prove to be a useful molecular mechanical tool for rational drug design. Especially in binding free energy calculations and molecular docking studies.

The Antechamber tool set is designed to allow the rapid generation of topology files for use with the AMBER simulation programs. This is useful in situations where you want to automatically screen a large number of compounds. Allowing antechamber to calculate charges and atom types automatically for us using GAFF allows it to be included in shell scripts that process a large number of compounds. However, as with any automated system it is not perfect and if you plan on focusing on a single system you should consider manually assigning atom types and carefully validating things. With Antechamber, one may solve the following problems:

1. Automatically identify bond and atom types

2. Judge atomic equivalence

3. Generate residue topology files

4. Find missing force field parameters and supply reasonable suggestions

Remember though that Antechamber is not a replacement for due diligence. You should always closely examine the atom types that Antechamber assigns and verify to yourself that the choices are reasonable. You should never use scientific software in a "Black Box" approach!

In this tutorial we shall use the Antechamber tools with LEaP to create topology and coordinate files for the prescription drug Sustiva (Efavirenz). Efavirenz) is a human immunodeficiency virus type 1 (HIV-1) specific, non-nucleoside, reverse transcriptase (RT) inhibitor marketed by Bristol Myers Squibb for controlling the progression of HIV infection in humans. The chemical name for Sustiva is  (S)-6-chloro-(cyclopropylethynyl)-1,4-dihydro-4-(trifluoromethyl)-2H-3,1-benzoxazin-2-one. Its empirical formula is C₁₄H₉ClF₃NO₂ and it's 2D structure is:

Here is a basic 3-dimensional geometry for Sustiva from which we will start to build our topology and coordinate files: sustiva.pdb. Our Sustiva pdb file is extracted from the RT-sustiva complex pdb file (PDB ID 1FKO ). The coordinates of sustiva are associated with a residue called "EFZ" (Efavirenz).
By all means open it up in VMD and take a look at it.

We shall use Antechamber to assign atom types to this molecule and also calculate a set of point charges. Antechamber is the most important program within the set of Antechamber tools. It can perform many file conversions and can also assign atomic charges and atom types. Depending on its inputs, antechamber executes the following programs (all provided with AmberTools): sqm, atomtype, am1bcc, bondtype, espgen, respgen and prepgen. It will also generate a series of intermediate files (all in capital letters).

Create parameter and coordinate files for Sustiva

First, let's run reduce to add all the hydrogen atoms to the pdb file. The hydrogenated sustiva coordinates can be found in sustiva_h.pdb.

reduce sustiva.pdb > sustiva_h.pdb  

(A note: reduce relies on a database of known ligands that are found in the PDB. If you have a molecule that is not already in the PDB component database, and you need to add hydrogens, give OpenBabel a try:

obabel sustiva.pdb -O sustiva_h.pdb -h  

Be sure to check that either reduce or obabel have made the correct changes. And note that OpenBabel is not distributed with Amber.)

To be consistent with the name of pdb, we will change the name of the residue from "EFZ" to "SUS",  and create a new pdb file: sustiva_new.pdb.

Now let's try using antechamber on our sustiva pdb file. To create the "mol2" file, required to define a new unit in LEaP, we simply run the following command:

antechamber -i sustiva_new.pdb -fi pdb -o sustiva.mol2 -fo mol2 -c bcc -s 2

Here the -i sustiva.pdb specifies the name of the 3D structure file and the -fi pdb tells antechamber that this is a pdb format file (we could easily have used any number of other supported formats including Gaussian Z-Matrix [gzmat], Gaussian Output [gout], MDL [mdl], amber Restart [rst], Sybyl Mol2 [mol2]). The -o sustiva.mol2 specifies the name of our output file and the -fo mol2 states that we want the output file to be of Tripos Mol2 format (this is a format supported by LEaP). The -c bcc option tells antechamber to use the AM1-BCC charge model in order to calculate the atomic point charges while the -s 2 option defines the verbosity of the status information provided by antechamber. In this case we have selected verbose output (2).

So, go ahead and run the above command. The screen output should be as follows:

Running: /usr/local/amber10/bin/bondtype -j full -i ANTECHAMBER_BOND_TYPE.AC0 -o ANTECHAMBER_BOND_TYPE.AC -f ac

Running: /usr/local/amber10/bin/atomtype -i ANTECHAMBER_AC.AC0 -o ANTECHAMBER_AC.AC -p gaff

Total number of electrons: 160; net charge: 0

Running: /usr/local/amber10/bin/sqm.sh

Running: /usr/local/amber10/bin/am1bcc -i ANTECHAMBER_AM1BCC_PRE.AC -o ANTECHAMBER_AM1BCC.AC -f ac
-p /usr/local/amber10/dat/antechamber/BCCPARM.DAT -s 2 -j 1

Running: /usr/local/amber10/bin/atomtype -f ac -p bcc -o ANTECHAMBER_AM1BCC.AC -i ANTECHAMBER_AM1BCC_PRE.AC

You should also get a whole series of files written to your directory.

ANTECHAMBER_AC.AC      ANTECHAMBER_AM1BCC_PRE.AC  ATOMTYPE.INF  sqm.out  sustiva.mol2
ANTECHAMBER_AC.AC0     ANTECHAMBER_BOND_TYPE.AC0  divcon.pdb    sqm.pdb  sustiva.pdb
ANTECHAMBER_AM1BCC.AC  ANTECHAMBER_BOND_TYPE.AC   sqm.in     

The files in CAPITALS are all intermediate files used by antechamber and are not required here. You can safely delete them. These files are not deleted by default since they may be of interest if things didn't work correctly. The sqm.xxx files are input and output from the sqm quantum mechanics code used by Antechamber to calculate the atomic point charges. We are not interested in the data here except to check that the sqm calculation completed successfully:

The file that we are really interested in, and the reason we ran Antechamber in the first place, is the sustiva.mol2 file. This contains the definition of our sustiva residue including all of the charges and atom types that we will load into LEaP to when creating our prmtop and rst7 files. Let's take a quick look at the file:

As you can see this file contains the 3 dimensional structure of our sustiva molecule as well as the charge on each atom, final column, the atom number (column 1), its name (column 2) and it's atom type (column 6). It also specifies the bonding at the end of the file. This file does not, however, contain any parameters. The GAFF parameters are all defined in $AMBERHOME/dat/leap/parm/gaff.dat. The other thing you should notice here is that all of the GAFF atom types are in lower case. This is the mechanism by which the GAFF force field is kept independent of the macromolecular AMBER force fields. All of the traditional AMBER force fields use uppercase atom types. In this way the GAFF and traditional force fields can be mixed in the same calculation.

While the most likely combinations of bond, angle and dihedral parameters are defined in the parameter file it is possible that our molecule might contain combinations of atom types for bonds, angles or dihedrals that have not been parameterised. If this is the case then we will have to specify any missing parameters before we can create our prmtop and rst7 files in LEaP.

We can use the utility parmchk2 to test if all the parameters we require are available.

parmchk2 -i sustiva.mol2 -f mol2 -o sustiva.frcmod

Run this command now and it will produce a file called sustiva.frcmod. This is a parameter file that can be loaded into LEaP in order to add missing parameters. Here it will contain all of the missing parameters. If it can antechamber will fill in these missing parameters by analogy to a similar parameter. You should check these parameters carefully before running a simulation. If antechamber can't empirically calculate a value or has no analogy it will either add a default value that it thinks is reasonable or alternatively insert a place holder (with zeros everywhere) and the comment "ATTN: needs revision". In this case you will have to manually parameterise this yourself. It is the hope that as GAFF is developed, the number of missing parameters will decrease. Let's look at our frcmod file:

remark goes here
MASS

BOND

ANGLE
ca-c3-c1   64.784     110.735   Calculated with empirical approach
c1-c1-cx   56.400     177.990   same as c1-c1-c3
c1-cx-hc   48.300     109.750   same as c1-c3-hc
c1-cx-cx   64.200     111.590   same as c1-c3-c3

DIHE

IMPROPER
ca-ca-ca-ha         1.1          180.0         2.0          General improper torsional angle (2 general atom types)
n -o -c -os        10.5          180.0         2.0          General improper torsional angle (2 general atom types)
c -ca-n -hn         1.1          180.0         2.0          General improper torsional angle (2 general atom types)
ca-ca-ca-n          1.1          180.0         2.0          Using default value

NONBON

We can see that there were a total of 4 missing angle parameters and 4 missing improper dihedrals. For the purposes of this tutorial we shall assume that the parameters Antechamber has suggested for us are acceptable. Ideally you should really test these parameters (by comparing to ab initio calculations for example) to ensure they are reasonable. If you see any parameters listed with the comment "ATTN: NEEDS REVISION" then it means that Antechamber could not determine suitable parameters and so you must manually provide these before you can proceed with the simulation. By default Antechamber will have set the values to zero.

We now have everything we need to load sustiva as a unit in LEaP. We just need to run tleap and ensure the GAFF force field is available.

$tleap -f oldff/leaprc.ff99SB

Once tleap is up and running we also need to ensure that it knows about the GAFF force field. There is a script in $AMBERHOME/dat/leap/cmd/ that will do this for us. We can load it into tleap with:

source leaprc.gaff

Our tleap console should now look something like this:

Now we can load our sustiva unit (sustiva.mol2): 

SUS = loadmol2 sustiva.mol2 

If you now type list in tleap you should see the new SUS unit (highlighted in bold): 

At this point we haven't loaded the frcmod file that parmchk2 gave us. Thus if we check our SUS unit we should find that there are 4 missing angle type parameters.

check SUS

Our missing angle type parameters were ca-c3-c1, c1-c1-cx, c1-cx-hc and c1-cx-cx. These correspond to the propyl ring and the c-c triple bond. This is what we would expect since this type of system is fairly rare in organic molecules. We can now load our frcmod file in order to tell tleap the parameters for these missing angle types.

loadamberparams sustiva.frcmod

If we now check out SUS unit we should find that there are no missing parameters:

We can now create the library file for sustiva (sus.lib), as well as the prmtop and rst7 files (sustiva.prmtop, sustiva.rst7).  

saveoff SUS sus.lib 

saveamberparm SUS sustiva.prmtop sustiva.rst7

The output from tleap shows a few warnings, which can be safely ignored (in this case!) due to the triangular bond geometry of sustiva: 

Instead of typing everything in the tleap console, the list of tleap commands mentioned above can also be used to create a tleap input file (tleap.in) and generate all the required files:

tleap -f tleap.in

Creating topology and coordinate files for Sustiva-RT complex

Since we can mix the traditional AMBER force fields with GAFF, we can at this point load a fragment of the reverse transcriptase (RT) from HIV virus and treat this using the ff99SB force field while treating the Sustiva molecule using the GAFF force field. We will need to use the Sustiva library file (sus.lib) that was created in the previous step.   

The RT-Sustiva complex can be found in the RCSB protein data bank (pdb code 1FKO). The corresponding pdb file is 1FKO.pdb.

HIV reverse transcriptase is a heterodimer composed of the p51 and p66 subunits. It is a large protein with a molecular mass of 117 kDa. For the purpose of this tutorial, we will use a truncated system in close proximity of Sustiva, including the finger and palm domains of the p66 subunit. The truncated pdb file is: 1FKO_trunc.pdb.

Now we are able to load the pdb file into tleap. 

First, we start tleap just like what we did in the previous section: 

tleap -f oldff/leaprc.ff99SB

>source leaprc.gaff 

>loadamberparams sustiva.frcmod

Now we load the Sustiva library file (sus.lib), followed by the complex pdb file 1FKO_trunc_sus.pdb.

>loadoff sus.lib 

>complex = loadpdb 1FKO_trunc_sus.pdb 

 

Again we can create a tleap input file (tleap2.in) and generate all the files for the Sustiva-RT complex:

tleap -f tleap2.in

Minimize and Equilibrate the Sustiva-RT complex

Once we have the topology and coordinate files of RT-Sustiva complex, we are ready to run short GB simulations on Sustiva-RT. Note, the procedure given here is very short in order to make the simulations compatible with the timescale of this tutorial. In a real "production" simulation you would typically run a much longer simulation (ns) in order to obtain good statistical convergence.

sander -O -i min.in -o 1FKO_sus_min.out -p 1FKO_sus.prmtop -c 1FKO_sus.rst7  -r 1FKO_sus_min.ncrst  &

You can use ambpdb to generate a pdb of the final minimized structures if you want:

ambpdb -p 1FKO_sus.prmtop -c 1FKO_sus_min.ncrst > 1FKO_sus_min.pdb

It is always a good idea to visualize your results since the human eye is very good at spotting anomalies.

The next step if to heat the RT-Sustiva complex. For speed we will do this very rapidly over 1ps. Ideally you should do this for much longer.

Here's our input file, we will run MD (imin=0) and this is not a restart (irest=0). In this example we will not use shake since it is possible that the hydrogen motion may effect the binding energy (probably not, but it serves as an example here). As we are not using shake we will need a time step smaller than normal. Here I will use a time step of 1 fs and run for 1000 steps [2 ps] (dt = 0.001, nstlim=1000, ntc=1). We will also write to our output file every 20 steps and to our trajectory [mdcrd] file every 20 steps (ntpr=20,ntwx=20). For temperature control we will use a Langevin dynamics approach with a collision frequency of 1 ps^-1. We will start our system at 0K and we want a target temperature of 300K (ntt=3, gamma_ln=1.0, tempi=0.0, temp0=300.0). And here's the input file: eq.in 

sander -O -i eq.in -o 1FKO_sus_eq.out -p 1FKO_sus.prmtop -c 1FKO_sus_min.ncrst  -r 1FKO_sus_eq.rst -x 1FKO_sus_eq.nc &

References

¹Wang, J., Wolf, R.M., Caldwell, J.W., Kollman, P.A., Case, D.A. "Development and Testing of a General Amber Force Field", J. Comp. Chem., 2004, 25, 1157 - 1173.

²Bayly, C.I., Cieplak, P., Cornell, W.D., Kollman, P.A. "A Well-Behaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges : The RESP Model", J. Phys. Chem, 1993, 10269-10280.

³Jakalian, A., Bush, B.L., Jack, B.D., Bayly, C.I., "Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model: I. Method.", J. Comp. Chem., 2000, 21, 132-146.