(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
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 waterwater 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 ionwater interactions would be included in the calculation and therefore alter the results for the watersolute 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 gistoutput.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 gistoutput.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 gistoutput.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 dTStransdens(kcal/mol/A^3) dTStransnorm(kcal/mol) dTSorientdens(kcal/mol/A^3) dTSorientnorm(kcal/mol) Eswdens(kcal/mol/A^3) Eswnorm(kcal/mol) Ewwdens(kcal/mol/A^3) Ewwnormunref(kcal/mol) Dipole_xdens(D/A^3) Dipole_ydens(D/A^3) Dipole_zdens(D/A^3) Dipoledens(D/A^3) neighbordens(1/A^3) neighbornorm ordernorm
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: dTStransdens(kcal/mol/Å^{3}), dTStransnorm(kcal/mol), dTSorientdens(kcal/mol/Å^{3}) and dTSorientnorm(kcal/mol).
c) Density weighted (dens) and normalized (norm) Solutewater and waterwater interaction energy: Eswdens(kcal/mol/Å^{3}), Eswnorm(kcal/mol), Ewwdens(kcal/mol/Å^{3}), Ewwnormunref(kcal/mol).
d) x, y, z dipole components, Dipole_xdens(D/Å^{3}), Dipole_ydens(D/Å^{3}), Dipole_zdens(D/Å^{3}), and the total dipole Dipoledens(D/Å^{3}).
e) Density weighted (dens) and normalized (norm) number of water neighbors within 3.5 Å of the voxel: neighbordens(1/Å^{3}) neighbornorm
f) Normalized order parameters, ordernorm.
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.
gistdipoledens.dx  Magnitude of mean dipole moment (polarization) (Debye/Å^{3}).
gistdipolexdens.dx  Xcomponent of the mean water dipole moment density (Debye/Å^{3}).
gistdipoleydens.dx Ycomponent of the mean water dipole moment density (Debye/Å^{3}).
gistdipolezdens.dx  Zcomponent of the mean water dipole moment density (Debye/Å^{3}).
gistdTSorientdens.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.
gistdTStransdens.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.
gistEswdens.dx  Density weighted SoluteWater interaction energy (kcal/mole/Å^{3}).
This is the interaction of the water with the entire solute. Both LennardJones 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.
gistEwwdens.dx  Density weighted WaterWater interaction energy (kcal/mole/Å^{3}).
This is the interaction of the water with the entire solvent. Both LennardJones 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.
gistgH.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.
gistgO.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.
gistneighbornorm.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.
gistordernorm.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 gistoutput.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 gistgO.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 bluewhitered scale. The following plot shows the isosurface of water density at the value of 3.0 gistgO.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 gistEwwdens.dx and gistEswdens.dx. These volumes represent volumes that contain either favorable solventsolute interactions or favorable solventsolvent interactions.
0.7 kcal/mol/Å^{3} isosurface of the solutewater interaction energy 
0.7 kcal/mol/Å^{3} isosurface of the waterwater 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, gistdTStransdens.dx and gistdTSorientdens.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 gistneighbornorm.dx and the 0.5 D isosurface of the total dipole gistdipoledens.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 gistgH.dx, water x, y, z dipole dipolexnorm.dx dipoleynorm.dx dipoleznorm.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 gistout.dat with the use of GistPP as shown in Section 4e.
(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