MBG wiki | RecentChanges | Blog | 2021-09-24 | 2021-09-23

Calculation of water and ion distributions

The object of the exercise is to generate a map showing the distribution of waters (hydration pattern) and ions around the solute. This map will pretty much look like a crystallographic electron density map, the difference being that instead of electron density it will show average water (or ion) density. An example of what we are trying to produce is shown below (taken from Reddy, Lecelrc & Karplus, Biophysical J, 84, 2003) :

 water and ions distribution around a DNA decamer

The contours depict the average water oxygen (yellow) and sodium (green) distributions around the MD-averaged structure of a DNA decamer.

Before starting with the practical issues, you should be warned about the two most notable shortcomings of the suggested procedure :

1. The method will fail if the solute diffused away from the water box. The reason for this is that the periodic boundary periodicity is not being taken into account (and the reason for this is programmer's lazyness).

The easiest way to check that the solute does stay inside the solvent box is to view the trajectory with carma (See also Quick view of large trajectory (DCD) files). For example, giving something in the spirit of

carma -gl 80.0 -trace -atmid " OH2 " -atmid " CA  " my.psf my.dcd

will pop-up a window showing the trajectory. Use the x, y, z keys to change the view. You should be seeing something like :

 Snapshots from carma showing a trajectory

2. The method will fail if the solute does not have a well behaved average structure. In other words, if the structure of the solute fluctuates too much, you should view the results with scepticism (or, repeat the calculation using shorter trajectory segments).

There are two ways of doing the calculation (both with carma). The first requires the creation of just one index file but creates a large intermediate DCD file. The second is doing the calculation on the fly (no intermediate files created) but requires you to create two index files (one for the atoms that will be used for fitting, the other for the atoms for which the density map will be build). You should probably use the second method.

Method 1 (creates large intermediate DCD file)

carma -v -fit -index -atmid " OH2 " -atmid " SOD " -atmid " CLA " -atmid " CA  " my.dcd my.psf

where the fit.index file (which should be present in the current directory) should contain indeces for the Cα atoms1.

structure @my.psf end
delete selection=( not (name OH2 or name CA or name CLA or name SOD)) end
write structure output=waters_ions.psf end
stop

 Orthogonal Coordinate limits in output file:
                 Minimum   Maximum    Centre     Range
       On X :     -24.30     24.30      0.00     48.60
       On Y :     -15.05     15.05      0.00     30.11
       On Z :     -10.98     15.87      2.44     26.85

carma -v -map -30 30 -20 20 -15 20 0.50 -atmid " OH2 " waters_ions.psf waters_ions.dcd
where the numbers after '-map' are the limits of the map and the grid size (in this case 0.50 Angstrom).

Method 2 (requires the creation of two index files)

structure @my.psf end
delete selection=( not (name OH2 or name CA)) end
write structure output=waters_CAs.psf end
stop

awk '{for (i=0 ; i < 112 ; i++) print i}'

 Orthogonal Coordinate limits in output file:
                 Minimum   Maximum    Centre     Range
       On X :     -24.30     24.30      0.00     48.60
       On Y :     -15.05     15.05      0.00     30.11
       On Z :     -10.98     15.87      2.44     26.85

carma -v -fit -index -fmap -30 30 -20 20 -15 20 0.50 -atmid " OH2 " -atmid " CA  " my.psf my.dcd

where the numbers after '-fmap' are the limits of the map and the grid size (in this case 0.50 Angstrom).

Viewing the map

The output from carma will be a map file named carma.fitted.dcd.map. This is an electron density format used by xplor and cns. For vmd to understand the format you will have to rename it so that its suffix is .edm. So :

mv carma.fitted.dcd.map carma.edm
vmd carma.edm

You will see nothing. Go to 'representations', select 'isosurface', select 'None' (instead of 'Box'), click 'create replication', and change the contouring level ('isolevel'). You should see it.

You can now quit vmd, use xplor to calculate an average solute structure (for the same frames) and then use vmd to view and study the model + map. If you have a crystal structure for your molecule, you could even compare the experimentally observed and the MD-derived hydration patterns.

 MD-derived hydration pattern example

 MD-derived hydration pattern of zwitterionic AFA tripeptide at 280K

One other interesting analysis you can perform is whether the distribution of waters (or ions) converged. A simple-minded way of doing this is to calculate the distributions for two non-overlapping parts of the trajectory (which will have to be a sufficiently long one). If the two maps (from the two trajectory segments) do look alike (say, have a high value of the linear correlation coefficient) then you can argue that the distribution converged.


Footnotes:

1. If you have, say, 112 Cα atoms, give

awk '{for (i=0 ; i < 112 ; i++) print i}'

and save the resulting lines in the fit.index file.