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
- 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.
- Protonate the file, e.g. by
protonate -k -d /home/bio/Data/Amber/PROTON_INFO.bio <pdb_file
>new_file_name
- Change the residue names, if needed (e.g. CYS to CY-; CYX to CYS,
HIP/HIE/HID to HIS)
- Remove crystal water molecules. It is normally best to do that,
except for metal sites.
- 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).
- cp /home/bio/Mead/mead-2.2.0p2/utilities/parse* .
(it copies to files, needed for the next step, to the current
directory).
- 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.
- 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.
- If you will use solinprot, move the solute to a separate file.
- 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
- -epsin is the dielectric constant within the molecule defined in
the file name.pqr.
- -blab2 gives intermediate much output. Alternatives are no
option (very little output), blab1, and blab3 (much output)
- -epssol gives the dielectric constant of the solution (i.e. the
rest of the space). Default = 80 = water.
- -ReactionField gives the reaction filed potential in a number of
points specified in the file name.fpt. The format for each
point
must be (x,y,z) (including the parenthesis). The units are such that
they
must be multiplied with a factor of 694.213 (econv*4.184/2) to
give
kJ/mole when multiplied with a charge in electron units.
- The result is in kcal/mole.
- The result is reproduced as the reaction field potential is
calculated in all atoms
Delta_G = econv/2 * Sum_over_all_atoms (reaction_field_potential *
charge)
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.
- You are not allowed to have a blank within the residue name
(therefore the "_" in the example above). In fact, the file is not read
with a fixed formate, but record-wise with blanks or tabs as delimiters.
- No atoms may have identical label and amino acid number -
There is no warning from the program if two atoms are identical, one of
them is just skipped.
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.
- The result depends slightly on the origin of the grid;
typically, the result can vary by 2.5 kJ/mole at different origins. It
can therefore be wise to do a few calculations where the grid origin is
displaced slightly and use the average as the best result. Then,
ON_ORIGIN is replaced by three real numbers separaded by blanks, e.g.
'0.1 0.1 0.1'.
- The grid should be extended beyond the protein. Gilson
recommends you to use the radius of the molecules*1.5. Gilson
recommends that a reasonable grid should be at least 1.5 times
as long as the protein in each
dimension (for solvation, 1.25 times may be enough), but our
calibrations
indicate that the result is not converged (to 1 kJ/mole) until you use
2.0 *the protein radius.
- Similarly, you seem to need a grid spacing of 0.1 Å
before the energy is converged to 1 kJ/mole. Gilson's recommendations:
For viewing field contours and lines, a grid spacin of 0.5-1.0
Å is usually adequate. For solvation energies, a grid sapcing
of about 0.2 Å is recommended. Bashford recommends a
grid spacing of ~0.25 Å for pKa calculations.
- The timing and memory demands of the program depend on the cube
of the number of grid points. Test what you can use on your mashine.
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)
- Calculate the potential of oxidised state in protein with
epsin=4, epsext=80 (ionicstr=0.1)
potential -epsin 4 -ionicstr 0.1 -blab2 name >name
.out
- 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.
- Use solvinprot to calculate the solvation energy of the ligand
in the protein with the metal.
The ligand is the solute, the protein+metal is the protein
solinprot -epsin1 4.0 -epsin2 4.0 -blab2 solute protein
- Calculate the solvation energy of the ligand in the
protein without the metal.
The ligand is the solute, the protein the protein
solinprot -epsin1 4.0 -epsin2 4.0 -blab2 solute protein
- 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