How to run MEAD 
(Macroscopic Electrostatics with Atomic Detail)

Mead is a program that is used for including solvation effects in biologcal systems such as proteins, etc. using an atomic model of the protein. This is done by solving the Poisson(-Boltzmann) equation in a dielectric medium, including distributed point charges, on a grid.
This is Poisson's equation (it relates the charge density, rho, and the electrostatic potential, phi. Epsilon is the dielectric constant, E is the electric field):
-grad * grad phi(r) = rho(r)/(epsilon*epsilon0) ( = div E(r) )
In the Poisson-Boltzmann equation, the right-hand side is supplemented with a term depending on mobile ions in the solvent. The equation will also be modified if the dielectric constant varies in the room (as it does in a protein); see the text by Gilson, below.

The current version is Mead-2.2.0 (1.1.8  and 1.1.5 are also available). 
The files are in /home/bio/Mead-2.2.0 and /home/bio/Mead/mead-2.2.0p2. 
To run the programs, add /home/bio/Mead/mead-2.2.0p2/Linux/bin (or Sgi/bin) and /home/bio/Mead/mead-2.2.0p2/utilities to your PATH
The program can be downloaded from http://www.scripps.edu/bashford .
Documentation (not very much): /home/bio/Mead/mead-2.2.0p2/README, /home/bio/Mead/mead-2.2.0p2/utilities/README, and /home/bio/Mead/Dok/dok.ps.
A good tutorial about continuum electrostatics and similar calculations (but with UHBD) can be found in http://gilsonlab.umbi.umd.edu (by Michael. K. Gilson; goto Intro Electrostatics).

Mead consists of several programs, three of which are described below (potential, solvate, and solinprot.

For solinprot, there is a detailed description by Torben Rasmussen.



How to start a simple calculation with Mead from a pdb file
  1. Start with a pdb file.
    If you have an equilibrated structure from Amber (with partial charges), changepdb command me may construct directly the pqr file. Then, you can continue from point 9 above after removing water molecules.  The Amber routine ambpdb, option -pqr also constructs a pqr file from Amber prmtop and coordinate files.
  2. Protonate the file, e.g. by
    protonate -k -d /home/bio/Data/Amber/PROTON_INFO.bio <pdb_file >new_file_name
  3. Change the residue names, if needed (e.g. CYS to CY-; CYX to CYS, HIP/HIE/HID to HIS)
  4. Remove crystal water molecules. It is normally best to do that, except for metal sites.
  5. With changepdb, set format to CNS (command f; 1; 3; 1; 1; 2; 1) and change the names of the HB3 atoms (command 3). Write out the file (command w).
  6. cp /home/bio/Mead/mead-2.2.0p2/utilities/parse* . 
    (it copies to files, needed for the next step, to the current directory).
  7. Run the Mead routine
    assign_parse_radch.pl <protonated_file >new_file.pqr 
    (note the prefix of the last file)
    It gives "PARSE" charges and radii of all atoms, except those listed by the program.
    Warnings of the type "Use of implicit split to @_ is deprecated at /assign_parse_radch.pl line 49." can safely be ignored.
  8. Manually, assign charges and radii of the atoms listed by the program. 
    Check that you have got integer charges of all residue with changepdb command cc.
  9. If you will use solinprot, move the solute to a separate file.
  10. Construct an ".ogm" file.
    Run changepdb, command xyz and take the three mid-points as the centre of the grid. Calculate the desired box size from the maximum value of diff. Decide what grid spacing you want to use (see below, ogm files) and calculate from this the number of grid points. Ensure that it is a odd number. You may use focusing, but only outside the protein.
If you have an equilibrate structure from Amber (with partial charges), changepdb command wpq may construct directly the pqr file. Then, you can continue from point 9 above.
The Amber routine ambpdb, option -pqr also constructs a pqr file from Amber prmtop and coordinate files.



Potential
You run it with:
potential -epsin 1 -epsext 80 -blab2 name > name.out
It calculates the potential of the molecule in the name.pqr file (using the grids, etc. in name.ogm) in a number of points given in the file name .fpt (format as above). The units is as above.
The result of solvate can be reproduced by two potential calculations, one with the same epsin and epsext=epssol as in solvate, the other with the same epsin but epsext=1. Then,
reaction_field = potential(epsin/epssol) - potential(epsin/1)

Files needed: name.ogm,  name.pqr, and name .fpt



Solvate
The program is executed with:

solvate -epsin 4.0 -blab2 name >> name .out

Files needed: name.ogm and name.pqr



Solinprot

solinprot -epsin1 4.0 -epsin2 4.0 -epsext 80 -blab2 solute protein >> solute.out

This program calculates the solvation energy of a molecule (solute; with epsin1) inside a protein (with dielectric constant epsin2), which together are emersed in a dielectric continuum of dielectric constant epsext. The solvation energy is calculated relative the molecule (with epsin1) in vacuum. The potential of protein point charges are included in the calculation.

Files needed: solute.ogm, solute.pqr, and protein .pqr
 

For solinprot, there is a detailed description by Torben Rasmussen.



The pqr file:
Generate the name.pqr file as a normal pdb file with charges and radii for all atoms after the pdb definitions. The format is actually free, and there must be a white space (blank) between each field. This can give errors for some protein atoms with four letters in some formats - be careful.
REMARK Test calculation of an imidazol; Merz-Kollman Charges obtained from QC calculation with Gaussian 98
REMARK SCF energy =    -3009.135870, it =   23, conv = 0.5931E-08, S**2 =  0.000
ATOM      8  N3  N_8     8      -1.940  -0.419  -0.594   -0.060685 1.50
ATOM      9  C4  C_9     9      -2.795   0.589  -0.560    0.072454 1.70
ATOM     10  N4  N10    10      -3.974   0.235  -1.138   -0.237111 1.50
ATOM     11  C5  C11    11      -3.866  -1.076  -1.564   -0.172007 1.70
ATOM     12  C6  C12    12      -2.597  -1.460  -1.216   -0.185532 1.70
ATOM     14  H1  N10    10      -4.786   0.829  -1.231    0.326699 1.00
ATOM     15  H2  C_9     9      -2.567   1.556  -0.123    0.112909 1.00
ATOM     23  H7  C12    12      -2.119  -2.418  -1.368    0.166500 1.00
ATOM     24  H8  C11    11      -4.677  -1.594  -2.054    0.172342 1.00
END
The pqr file can be created from an Amber-pdb file with charges using
changepdb, option mead.
This utility adds PARSE charges to each atoms [Sitkoff, Sharp and Honig, 1994, J. Phys. Chem. 98:1978--1988]. In practise, the radii a constructed to reproduce the data base file: /home/bio/Mead/MEAD-1.1.8/utilities/parse_radii .
There is a pear-script in the same directory assign_parse_radch.pl that is said to add both radii and charges to a pdb file, but I have not managed to run it properly.

Mats has done a program
/home/bio/Mats/Gausstopqr/gausstopqr
which convert a Gaussian MK calculation to a pqr file. I have not tested it.



The ogm file
The name.ogm file contains information of the grid. The first item is the definition of the origin (it could be useful to have it "on the origin", the second item is the number of grid lines (it has to be an odd number), and finally the spacing between the grid lines. Normally, this file should contain only one line, but sometimes you can focus the calculation on the system of interest with a series of values, the coarser ones only giving proper boundary conditions to the finer ones. Use this with great care, however. Example:
ON_ORIGIN           27    1.0
ON_ORIGIN           51    0.5
ON_ORIGIN           99    0.25
ON_ORIGIN          193    0.125


The fpt file
It contains the points for which the potential should be calculated. They are given as a row-list of coordinates in the format
(0.0, 0.0, 0.0)
(2.0, 0.0, 0.0)

where the parenthesis and commas need to be written out explicitly.

changepdb, command fpt constructs a .fpt file from a pdb (or pqr file) (all atoms included).


Test
You can reproduce exactly the Born formula using one atom:
ATOM      1  CU  CU1     1       0.000   0.000   0.000    2.000000 2.00
and
solvate -epsin 1 name
which gives -328.5 kcal/mole, compared to Born (q=2e, r=2 Å) -327.9 kcal/mole.



To compile on SGI
Add
DEFS=-DUSE_RTTI
in Makefile.sgi_IRIX6



To do a reorganisation energy calculation

1. Simplest model (ignoring the electronic polarisation and conformation change)

  1. Calculate the potential of oxidised state in protein with epsin=4, epsext=80 (ionicstr=0.1)

  2. potential -epsin 4 -ionicstr 0.1 -blab2 name >name .out
     
  3. Calculate the reorganisation energy by multiplying the resulting potentials with the charges of the oxidised state minus those of the reduced site and sum over all atoms in the site. This can be done with changepdb, option POtential.


To do a relative solvation energy calculation
e.g. how does a metal affect the binding of another ligand.
  1. Use solvinprot to calculate the solvation energy of the ligand in the protein with the metal.

  2. The ligand is the solute, the protein+metal is the protein
    solinprot -epsin1 4.0 -epsin2 4.0 -blab2 solute protein
     
  3. Calculate the solvation energy of the ligand in the protein without the metal.

  4. The ligand is the solute, the protein the protein
    solinprot -epsin1 4.0 -epsin2 4.0 -blab2 solute protein
     
  5. The desired energy is the difference of these two energies


Subject: CCL:MEAD pKa calculations
Date: Thu, 13 Feb 2003 12:32:32 -0800
From: Don Bashford <bashford@scripps.edu>
 

>>>>> "CUI" == CUI, Guanglei <cuigl@morita.chem.sunysb.edu> writes:

    CUI> Dear CCL readers, I'm trying to use MEAD to calculate pKa of
    CUI> three residues in a ligand binding site. The structure has
    CUI> been minimized with AMBER force field with explicit
    CUI> solvent. The charges of these residues have been derived from
    CUI> RESP calculations. All needed files are created and multiflex
    CUI> calculation used epsin = 4.0. After the calculation, the
    CUI> intrinsic pKas of two residues are very surprising. So here's
    CUI> my questions:

Hopefully you've used a fine grid of 0.25 Ang or so centered on
"center of interest".

    CUI>      1. What could be reasons that the self interactions for
    CUI> residue UNK and LYS are very big.

Probably because the accessibility of these residues to solvent is
low.  This tends to disfavor the charged states, and therefore lowers
pKint for Lys, and raises it for things like Asp and Glu.


    CUI> 2. During the multiflex
    CUI> calculation, I was warned of "vertex found with count = 2"
    CUI> three times. How should I deal with this?

This is almost certainly harmless.

    CUI> 3. Do I need to
    CUI> include explicit solvent molecules in the calculation,
    CUI> generally speaking? The binding site is not completely
    CUI> buried.

People have done this both ways.  It can be argued that keeping
crystallographic waters is more "detailed", but on the other hand,
when the ionization state of a residue changes these waters may move
or reorient, and a calculation that regards them as having fixed
coordinates will not capture this.  A blob of high dielectric can at
least repolarize!  In my experience, ordinary protein residue pKa's
are not much improved by including crystallographic waters, and may
even get a bit worse.  However, for multivalent things like metal
centers, an explicit hydration shell helps, because the hydration
shell is tightly bound and a singe change change doesn't really
re-orient it so much.

    CUI> 4. LYS (protonated) has 3 equivalent hydrogens. It's
    CUI> not very clear to me which one should be zero charged in
    CUI> LYS.st file.

My guess is it won't make so much difference, but why not try it all
three ways?

    CUI>      Here're the .summ file. Sorry for the lengthy
    CUI> description. Any suggestions are welcome. Thanks in advance.

    CUI>  site name           pKmod      delta self    delta back      pkint
    CUI>    UNK-259           7.8       23.1134      -5.41211       25.5013
    CUI>    LYS-162          10.4      -11.4993       1.32968        0.23042
    CUI>    TYR-155           9.6      0.694971      -1.28149        9.01348

What is UNK?  I'll assume it's negatively charged.  It's delta self is
indeed quite big.  It may be deeply buried.  It's charge may be very
compact or large (-2?).  The delta self you get is rather typical of a
buried lysine.  Do UNK and Lys interact strongly?  (see the .g file)
If so, the favorable charge-charge interaction may stabilize the
charge pair in spite of the large delta self terms.  I've seen such
things happen, for example in bR (Spassov et al., 2001, JMB 312:203).

Cheers,
Don Bashford




Mats 0006, Ulf 0008, 0105