(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Ross Walker 2015

AMBER Advanced Tutorial 25
SECTION 3

GIST Analysis Example
Analysis of water thermodynamics using Grid Inhomogeneous Solvation Theory of Factor Xa active site.

By Romelia Salomon, Crystal Nguyen, Steven Ramsey, Jonathan D. Gough, Mike Gilson, Tom Kurtzman & Ross Walker

3) Analyzing the MD results with GIST

The final stage of our calculation is to use the cpptraj implementation of GIST to calculate the solvent structural and thermodynamic properties of our system using the results of our MD simulation.

As the name implies (Grid Inhomogeneous Solvation Theory), a grid is used to define the region about which the solvent properties are calculated. Therefore the first step is to define the region of interest, in space via coordinates.

The following images show different views of the grid box employed for our GIST calculation. In order to understand the water structure about the ligand binding site, it is necessary that the box is placed exactly on top of the region where the ligand would bind. This is important as it defines the space where the thermodynamic quantities will be estimated. One option is to have a grid almost as large as the periodic box, however the user should be aware that this choice will make the GIST calculations much slower.

The images below show the grid as defined in this tutorial for Factor Xa active site in different orientations. The first figure shows the grid on the side with the ligand inside for visual aid. The last figure shows the box viewing the active site from the front as reference for the results that will be presented next.

Grid in side orientation plus ligand

Grid in side orientation without ligand

Grid in front orientation without ligand


 

For this tutorial it is assumed that you have downloaded cpptraj and know the basics of how to use it. The command line for GIST is simple and requires few parameters.

The cpptraj command options are as follows.

gist [doorder] [doeij] [refdens <rdval>] [gridcntr <xval> <yval> <zval>] [griddim <xval> <yval> <zval>] [gridspacn<spaceval>] [out <filename>]



The most important parameters determine the region being considered for analysis and therefore determine the amount of computational time required. They are as follows:

gridcntr: Coordinates (Å) of the center of the grid. The default values are (0.0,0.0,0.0).
griddim: Integer number of grid increments along each coordinate axis. The default values are (40,40,40).
gridspacn: Grid spacing (linear dimension of each voxel) in Angstroms. Values greater than 0.75 Å are not recommended, the default value is 0.5 Å .



The other optional paramaters provide additional thermodynamic and structural information, they are as follows:

doeij: Specifies if the triangular matrix representing the water-water interactions between pairs of voxels is calculated (default is false, not calculated).
dorder: Specifies whether the water order parameter will be calculated for each voxel (default is false, not calculated).
redens: Specifies a reference density of bulk water, used in computing g_O, g_H, and the translational entropy. The default value is 0.0334 molecules/Å3.



The present implementation of GIST supports the following water models:

a) TIP3P
b) TIP4P
c) TIP4PEW
d) TIP5P
e) TIP3PFW
f) SPCE
g) SPCFW

You should read the AmberTools manual to see the full set of parameters and available options for GIST.



For the analysis of the Factor Xa active site, we remove the ions prior to invoking GIST as the ion-water interactions would be included in the calculation and therefore alter the results for the water-solute interaction energy. The trajectory is large and therefore the GIST calculation takes several hours.

For this analysis we use the following cpptraj commands:

parm ChA_IONS_disul_TIP3P.prmtop
trajin prod_100ns_aligned.nc
strip @Ca
strip @Cl-
gist gridcntr 25.0 31.0 30.0 griddim 41 41 45 gridspacn 0.50 out gist-output.dat
go
quit

Simply, this script (gist.ptraj) can be run as follows:

cpptraj -i gist.ptraj > gist.ptraj.out &

As you will notice, a grid spacing value of 0.50 Å was employed. We do not recommend using values larger than 0.75 Å. Although, larger values will provide smoother results, the estimation of the thermodynamic values will be inaccurate. GIST is designed to have at most one water molecule per voxel for each frame.

In our example, GIST produces the file gist-output.dat, which provides all calculated quantities in a tabular format, where each row corresponds to one grid box (voxel).

Below, we show the first three lines: the default comment, the column headers, and the first line of numerical output. The gist-output.dat file is a collection of all data produced by GIST printed in the order of grid point.

GIST Output, information printed per voxel

voxel xcoord ycoord zcoord population g_O g_H dTStrans-dens(kcal/mol/A^3) dTStrans-norm(kcal/mol) dTSorient-dens(kcal/mol/A^3) dTSorient-norm(kcal/mol) Esw-dens(kcal/mol/A^3) Esw-norm(kcal/mol) Eww-dens(kcal/mol/A^3) Eww-norm-unref(kcal/mol) Dipole_x-dens(D/A^3) Dipole_y-dens(D/A^3) Dipole_z-dens(D/A^3) Dipole-dens(D/A^3) neighbor-dens(1/A^3) neighbor-norm order-norm

0 15 21 19 415 0.994012 0.947305 0.000118875 0.00358058 0.000822792 0.0247829 0.000326025 0.00982002 -0.318998 -9.60836 0.00262741 -0.00135397 0.00170226 0.00341089 0.18016 5.42651 0

The first entries correspond to the voxel number followed by its x, y, z coordinates, while the next entries specify different values:

a) Population, atom density for oxygen and hydrogen (g_O and g_H).
b) Density weighted (dens) and normalized (norm) Translational and Orientational entropy: dTStrans-dens(kcal/mol/Å3), dTStrans-norm(kcal/mol), dTSorient-dens(kcal/mol/Å3) and dTSorient-norm(kcal/mol).
c) Density weighted (dens) and normalized (norm) Solute-water and water-water interaction energy: Esw-dens(kcal/mol/Å3), Esw-norm(kcal/mol), Eww-dens(kcal/mol/Å3), Eww-norm-unref(kcal/mol).
d) x, y, z dipole components, Dipole_x-dens(D/Å3), Dipole_y-dens(D/Å3), Dipole_z-dens(D/Å3), and the total dipole Dipole-dens(D/Å3).
e) Density weighted (dens) and normalized (norm) number of water neighbors within 3.5 Å of the voxel: neighbor-dens(1/Å3) neighbor-norm
f) Normalized order parameters, order-norm.



GIST produces a series of dx files in addition to the tabular output file, which contains all of the calculated data. These files can be readily visualized using VMD.

gist-dipole-dens.dx - Magnitude of mean dipole moment (polarization) (Debye/Å3).

gist-dipolex-dens.dx - X-component of the mean water dipole moment density (Debye/Å3).

gist-dipoley-dens.dx- Y-component of the mean water dipole moment density (Debye/Å3).

gist-dipolez-dens.dx - Z-component of the mean water dipole moment density (Debye/Å3).

gist-dTSorient-dens.dx - Density weighted first order orientational entropy (kcal/mole/Å3).
The more negative the isovalue, the more restricted or unfavorable the orientation of the water is. The value at bulk density is expected to be zero, most values will fall between 0 and -1.
gist-dTStrans-dens.dx - Density weighted first order tranlational entropy (kcal/mole/Å3).
The more negative isovalue, the more restricted the water is positionally. The value at bulk density is expected to be zero, most values will fall between 0 and -1.
gist-Esw-dens.dx - Density weighted Solute-Water interaction energy (kcal/mole/Å3).
This is the interaction of the water with the entire solute. Both Lennard-Jones and electrostatic interactions are computed without any cutoff, within the minimum image convention. The more negative the number, the more favorable the interaction between the water and solute is. Isovalues will be negative numbers.
gist-Eww-dens.dx - Density weighted Water-Water interaction energy (kcal/mole/Å3).
This is the interaction of the water with the entire solvent. Both Lennard-Jones and electrostatic interactions are computed without any cutoff, within the minimum image convention. The more negative the number, the more favorable the interaction between water within this voxel and all other water molucules is. Isovalues will be negative numbers.
gist-gH.dx - Number density of hydrogen centers found in the voxel, in units of the bulk density.
The expectation value of g_H for a neat water system is unity. The units of this dx file are the density/bulk density. Therefore an isovalue of 1 represents every voxel at or greater than bulk density.
gist-gO.dx - Number density of oxygen centers found in the voxel, in units of the bulk density.
The expectation value of g_O for a neat water system is unity. The units of this dx file are the density/bulk density. Therefore an isovalue of 1 represents every voxel at or greater than bulk density.
gist-neighbor-norm.dx - Mean number of neighboring water molecules, per water molecule found in the voxel (units of number per water).
Two waters are considered neighbors if their oxygens are within 3.5 Angstroms of each other.
gist-order-norm.dx - Average Tetrahedral Order Parameter for water molecules found in the voxel, normalized by the number of waters in the voxel.
doorder must be declared in the command line for this file to be produced.

All the information GIST produces allows for the quantitative thermodynamic analysis of each point of the grid. Information is produced on population and normalized density which show the level of water occupancy of that position (Note: normalized values are not written to dx files, they are available in the gist-output.dat file. For more information see Section 4e in the next part of the tutorial). Values such as interaction energy and entropy show if a water molecule being present at that position is favorable or not with respect to bulk water conditions. This information is valuable for understanding solvation and binding.

VMD can be used to visualize the dx files on the protein surface. To do so:

Load VMD
Load your protein structure (ChA_IONS_disul_TIP3P.prmtop followed by the ChA_IONS_disul_TIP3P.inpcrd choosing Amber7 Restart as the file type
In the "Graphical Representations" panel, enter "all not water" ino the "selected Atoms" field
In the "Graphical Representations" panel: in the "Coloring Method" select "Color ID" and change the value to #2 (Grey)
In the "Graphical Representations" panel: in the "Drawing Method" select Surf."
You should now see the protein as a grey surface

Select Load files for: New molecule
Click Browse and select the gist-gO.dx file
Click Load
you should now see a grid box with voxel points

In the "Graphical Representations" panel: in the "Draw" menu, select "Wireframe"
In the "Graphical Representations" panel: in the "Show" menu, select "Isosurface"
In the "Graphical Representations" panel: in the "Isovalue" menu, enter in 3.0 to visualize voxels at 3 x bulk density
you should now see an image similar to the first figure below

These instructions can be utilized to produce a visualization similar to the figure below, but can also be applied to any dx file output from Gist using isovalues appropriate for that .dx file (see above, or reference the figures presented below for examples).


Here we show some plots produced with GIST output at a given isosurface on top of the electrostatic potential energy surface of the protein, in a blue-white-red scale. The following plot shows the isosurface of water density at the value of 3.0 gist-gO.dx. This figure highlights volumes where water molecules are located during most of the simulation.

The following plots represent isovalues of -0.7 (kcal/mol/Å3) for both gist-Eww-dens.dx and gist-Esw-dens.dx. These volumes represent volumes that contain either favorable solvent-solute interactions or favorable solvent-solvent interactions.

-0.7 kcal/mol/Å3 isosurface of the solute-water interaction energy

-0.7 kcal/mol/Å3 isosurface of the water-water interaction energy



Similar plots can be generated to represent regions with unfavorable entropy. The following plots show -0.1 (kcal/mol/Å3) isosurface's of density weighted orientational and translational entropy, gist-dTStrans-dens.dx and gist-dTSorient-dens.dx.

-0.1 kcal/mol/Å3 isosurface of orientational entropy

-0.1 kcal/mol/Å3 isosurface of translational entropy



The following plot shows the 5.6 isosurface of the neighbor gist-neighbor-norm.dx and the 0.5 D isosurface of the total dipole gist-dipole-dens.dx.

5.6 isosurface of the neighbor

0.5 D/Å3 isosurface of the total dipole

Here we include the rest of the files produced by GIST: Hydrogen total density gist-gH.dx, water x, y, z dipole dipolex-norm.dx dipoley-norm.dx dipolez-norm.dx, respectively. If you are feeling adventurous you could now try repeating this with a longer simulation, like the one used in the paper [Nguyen et al.].

One can extract additional GIST data from gist-out.dat with the use of GistPP as shown in Section 4e.


CLICK HERE TO GO TO SECTION 4


(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Ross Walker 2015