> ---- Karen Callahan wrote ----:
> The easiest thing to do is look up the OPLSAA force
> field, it is pretty good for Mg2+, so it might be good
> for Ca2+.
The vdw radius for Ca2+ in OPLSAA is 1.3537 (derived from its sigma value; no indication of source).
> You may also try the Charm force field for r-Ca2+.
It is 1.3670 in CHARMM (from the sigma value in Marchand & Roux 1998 Proteins 33: 265-284, which is derived to reproduce experimental Ca2+ hydration free energy).
These two values are close to Aqvist's 1.3263 (Aqvist 1990 JPC 94: 8021-8024), which is derived to reproduce experimental Ca2+ hydration free energy and Ca2+ to water oxygen distance in SPC water (results transferrable to that in TIP3P water). The parm99.dat file indicates that Ca2+ vdw parameters (atom type C0) are from Aqvist's work, however, it is 1.7131 in parm99.dat. Can anyone tell me why?
> You will want to come up with some way to verify that
> whichever force field you use is giving reasonable
> results.
The Aqvist values give much more reasonable Ca2+ coordination than the default parm99.dat values in my implicit-solvent minimization of the crystal structure (see my previous email below). However, I also notice that the average distance between Ca2+ and 6 coordinated oxygens is reduced from experimental 2.40 to 2.30 after a 400-SD minimization with the Aqvist values. The underestimation of calcium-oxygen distance has also been observed in the Marchand & Roux's work above. They claim it is due to the lack of nonadditive polarization in partial charges. I also agrees with Bill Ross's point that the Ca2+ vdw radius is underestimated during its derivation in TIP3P/SPC water models that have an expanded vdw for oxygen and no vdw for H (
http://amber.ch.ic.ac.uk/archive/200504/0377.html ). The water oxygen's vdw radius is about 0.1 larger than those protein oxygens (atom type O and O2) in my system, so I added 0.1 to Aqvist's Ca2+ vdw radius without chaing the epsilon value and re!
-m!
inimized the system. The minimized structure has an ave
2.38 between Ca2+ and the 6 oxygens, which is close to the experimental structure's 2.40 and idential to the statistical value surveyed on Cambridge Structural Database by Harding (Harding 1999 Acta Cryst. 1999 D55: 1432-1443). But, without doubt, this arbitary value will not reproduce Ca2+'s interaction with TIP3 waters anymore.
Therefore, I think the best way to derive the vdw parameters for structural and solvent-accessable Ca2+ is to saitsfy the following conditions:
1. Reproducing experimental Ca2+ hydration free energy
2. Reproducing experimental average distance between Ca2+ and water oxygens
3. Reproducing experimental interaction energy between Ca2+ and protein oxygens (but I doubt such experimental data exists)
4. Reproducing experimental average distance between Ca2+ and protein oxygens
Aqvist's values have satisfied condition #1 and #2, but as Bill Ross has pointed out that there is a range of vdw values that will yield the same radial distribution and solvation free energies in TIP3 water (
http://amber.ch.ic.ac.uk/archive/200504/0377.html ). Thus, one can use conditions #3 and #4 to specify the vdw values. Has anyone done this? Is it worthy to do that or one can just live with the ~0.1 underestimation of the average distance between Ca2+ and protein oxygens since there are other inaccuracies in the force field?
Another concern I have is that Ca2+ is not the only atom type (C0 in parm99.dat) whose vdw parameters are derived in SPC/TIP3P water models. For example, the 6 oxygens in my systems have atom types O and O2, whose vdw parameters are from the OPLS force field and derived in SPC/TIP3P water models as well (see OPLS development in Jorgensen & Tirado-Rives 1988 JACS 110: 1657-1666). So again, is it worthy to adjust Ca2+ vdw parameters alone?
Thanks and happy holiday.
Sincerely,
Qing
----- Previous Message ----
From: Qing Zhang <qingzhang_nyu.yahoo.com>
To: amber.scripps.edu
Sent: Friday, December 8, 2006 4:31:02 PM
Subject: Re: AMBER: Simulating proteins with calcium ion
Thanks
for the r
positions, and its points to the van der Waals parameters
for Ca2+ in parm99.dat. I am giving a more detailed description of the
system/problem and my analysis. The description will automatically
answer some questions by Fenghui and Tom, and I will explicitly anwer
the rest.
The system:
A protein-protein complex with a
structural Ca2+. The Ca2+ ion binds to 6 oxygens from 5 residues (2
GLN, 2 ASP, and 1 GLY). The 6 oxygens form a binding pocket like a
half-sphere, and Ca2+ is located nearly at the sphere center. The
distances between the oxygens and Ca2+ range from 2.2 to 2.5 Angstrom.
The problem:
Energy
minimizations of the complex crystal structure (using AMBER 9) cause
Ca2+ to escape the half-sphere by about 1.9 Angstrom to the solvent but
still binds to a few oxygens (Tom, there is no move/penetration through
other atoms). I used 400 SD followed by 600 CG with igb=5 and a 3
kcal/mol (or 10) restraint on the heavy atoms of the proteins and Ca2+.
Crystal waters are built into the system but there is no explicit
sovlent. A sample input file is below.
--------------------------------------------------------
Minimization with Cartesian restraints
&cntrl
imin=1, maxcyc=1000, ntmin=1, ncyc=400,
scee=1.2,
igb=5,
ntb=0, cut=16.0,
ntpr=100,
ntr=1, restraint_wt=3.0, restraintmask=':1-297 & !.H=',
ntc=2,
/
--------------------------------------------------------
The analysis:
The
first thing came to my mind is the radius of Ca2+ in parm99.dat. As the
6 oxygens (type O2 and O) have radii of 1.6612 (parm99.dat) and the
distances between Ca2+ and the oxygens range from 2.2 to 2.5, the large
raius for Ca2+ in parm99.dat (1.7131) will cause a large van der Waals
penalty and make Ca2+ to be pushed out of the binding pocket. So I did
energy composition analysis on Ca2+ (idecomp=2 and the restraint is
removed to make idecomp work). The initial structure with only 1 SD
gives a vdw of
nal |vdw |eel |pol |sas
290 0.000 49.029 -353.216 32.
57 0.000
Minimization
with 400 SD (no CG as Tom suggested it might cause large jumps) leads
to the same dislocation of Ca2+ and gives a vdw of 6 and eel of -350:
resid |internal |vdw |eel |pol |sas
290 0.000 5.994 -350.108 19.389 0.000
It indicates that Ca2+ is pushed out of the binding pocket during minimization to reduce van der Waals penalty.
In
order to further prove it, I reduced the radius of Ca2+ to about 1.3.
This is based on a post by Kenley Barrett on AMBER achieve. In this
post, the vdw parameters for divalent ions computed based on the Aqvist
paper (JPC 1990,k 94: 8021) are listed:
http://amber.ch.ic.ac.uk/archive/200504/0376.html
I
took the vdw parameters of Ca2+ (1.3263, 0.4497) from the post,
replaced those (for C0) in parm99.dat, and re-generated the topology
file. The energy decompostion on the initial structure with 1 SD show a
vdw of only 6 (reduced from 49):
resid |internal |vdw |eel |pol |sas
290 0.000 6.452 -353.216 32.657 0.000
Then
I minimized the system for 400 SD (same minimization condition as the
run without radius modification). The Ca2+ ion stays at its crystal
location! The movement is only 0.13 Angstrom, and the energy
decomposition shows a vdw of 12 and eel of -396:
resid |internal |vdw |eel |pol |sas
290 0.000 11.906 -395.934 61.374 0.000
From
these observations, the instability of Ca2+ during the minimizations is
clearly due to the vdw parameters for Ca2+ in parm99.dat.
I am not familiar with divalent vdw parameterization. If someone has more insights on this, please feel free to raise them.
Thanks,
Qing
======================================
Qing Zhang, Ph.D.
Research Associate
Department of Molecular Biology, MB-5
The Scripps Research Institute
10550 North Torrey Pines Road
La Jolla, CA 92037-1000
Tel: (858) 784-2333
Fax: (858) 784-2860
ebsite: http://www.qingzhang.info
======================================
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Sun Dec 24 2006 - 06:07:41 PST
