[gmx-users] pKa calculation from MD trajectories

[baptista at itqb.unl.pt]

Dear Peter, David and Ivano,

First, I note that running an MD simulation and doing pKa calculations
*afterwards* (using the trajectory conformations) yields something which is
not well defined from a statistical thermodynamic viewpoint. There are
several ad hoc methods that yield "conformation-averaged" pKa's in this
way, but in all cases the conformation and protonation states are
incorrectly coupled, as I discussed elsewhere [Baptista et al (2002)
J. Chem. Phys. 117:4184-4200].

However, if the pKa calculations are done *during* the MD run *and* used to
feed it with new charges, then one may get what is technically known as
constant-pH MD. The aim of constant-pH MD is not merely to mix MD and pKa
calculations, but rather to do simulations that jointly sample protein
conformations *and* protonation states in a correct way for the given pH
("correct" meaning "from the proper statistical thermodynamic ensemble").
Thus, the solution pH becomes a simulation parameter like the temperature
and pressure. The idea was first put forward and applied to BPTI in a paper
I wrote some years ago [Baptista et al (1997) Proteins 27:523-544],
actually using TITRA; maybe this is one of the papers Peter was referring
to.

When performing constant-pH MD one should be aware of three important
points: (1) how are protonation energies computed; (2) how are protonation
states sampled; (3) how is the protonation data inserted into the MD run.

Point (1) requires solving the Poisson-Boltzmann (PB) equation for the
protein and the surrounding solvent. A simple approach is to assume a
spherical protein and use the modified Tanford-Kirkwood (MTK) model, as
implemented in TITRA by Paulo Martel. But, although this model is very
fast, the neglect of the protein shape and of the self-energy terms may has
serious consequences, as was clearly shown many times by Warshel, Honig,
and others. The solution is to use a continuum electrostatic (CE) model
with atomic detail, and solve the PB equation numerically, as implemented
in DELPHI, UHBD, MEAD, etc (an alternative is Warshel's PDLD model).
Although I used TITRA some years ago for performance reasons, I think that,
now that computer power has increased significantly, we can all move on to
these more realistic methods.

With respect to point (2), enumerating all protonation states is usually
unfeasible and one must resort to approximations. The simplest one is to
use a mean field (MF) approximation, ie, consider that the sites see an
average state of each other; this is the approach of the MTK model
implemented in TITRA.  However, since the MF approximation leads to
incorrect predictions when strongly interacting sites titrate in the same
pH region, it is safer to use less restrictive methods, such as partial MF
approximations or Monte Carlo (MC) sampling. I tend to use MC, though one
should always be alert to sampling problems.

As for point (3), I have proposed two approaches: the implicit titration
(IT) method [Baptista et al (1997) Proteins 27:523-544], and the stochastic
titration (ST) method [Baptista et al (2002) J. Chem. Phys. 117:4184-4200].
In IT method the titrable sites have average charges that are updated
during the MD simulation, while in ST the sites change their (non-average)
charges in a stochastic way. The ST has some theoretical and practical
advantages over IT, discussed in that paper. The application of ST to
succinic acid worked quite well, and more work is in progress for peptides
and proteins. The current implementation consists of a slightly modified
version of GROMOS96, MEAD, and my own MC program (MCRP). The important
thing is that the MD performance is not drastically affected. Implementing
it in GROMACS is envisioned, given its speed advantages.

Finally, I should note that two other methods were recently proposed for
constant-pH MD, *not* using CE [Borjesson & Hunenberger (2001) J. Chem.
Phys. 114:9706-9719; Burgi et al (2002) Proteins 47:469-480]. Although I
have some doubts on the validity of the first [Baptista (2002) J. Chem.
Phys.  116:7766-7768], they show a renewed interest in constant-pH MD, and
I am confident that, in a form or another, it will become a standard
simulation tool in the near future.

Best regards,
Antonio

-- 
Antonio M. Baptista
Instituto de Tecnologia Quimica e Biologica, Universidade Nova de Lisboa
Av. da Republica, EAN, ITQB II, Piso 6, Apartado 127
2781-901 Oeiras, Portugal
phone: +351-214469613         email: baptista at itqb.unl.pt
fax:   +351-214411277         WWW:   http://www.itqb.unl.pt/~baptista
--------------------------------------------------------------------------

>On Mon, 2003-04-28 at 12:53, Peter Fojan wrote:
Dear David,
the program is tested on a broad range of proteins and there are a number of
papers available on the electrostatic calculations on proteins including on
where this approach has been used to monitor the pk related changes on a peptide
during an MD run. The program performs actually very nice, I dont know about
MEAD but it has been tested against UHBD. The output from this program can be
used a input files for DELPHI and the charge distribution can be visualized
in GRASP.
The program is available from us for non profit organisations free of any charges.
The only thing to do is to fax us a signed license agreement and I deliver the
executables by email.
The only restriction so far is that it can only run under SGI IRIX. I still
have to port it to standard linux for a broader availability. It also has the
possibilities of taking measrued pKa values from NMR experiments and include
these in the calculation.
I hope this piece of information is useful,
Best,
    Peter
PS: I am working on a LINUX port of the program in my scare spare time. When
there is a release ready I will elt everybody who is interested know.

>> Dear Ivano,
>> we have developed a program called TITRA that can calculate the Pka
>> values of titratble residues in a protein. The input files would be the
>> structure of the protein. So basically this one could do it, but it
>> would be a bit tedious, becasue you would need snapshots from your
>> trajectory and process every single one with titra. But basically it can

>> be done.
>
>And it should! That would be worth a paper I'd say. Is your program
>available somewhere? How does it compare to MEAD?
>
>> Best,
>>       Peter
>>
>> Ivano Eberini wrote:
>>
>> > Hi,
>> > we run some MD simulations for beta-lactoglobulin, starting from
>> > slightly different conditions.
>> > We were wondering if it is possible to calculate the pKa value for a
>> > specific residue (a glutamic acid) from MD trajectories.
>> > Thanks in advance for any suggestion and best regards,
>> >
>> > Ivano Eberini
>> > --
>> > Ivano Eberini
>> > Gruppo di Studio per la Proteomica e la Struttura delle Proteine
>> > Dipartimento di Scienze Farmacologiche
>> > Università degli Studi di Milano
>> > Via G. Balzaretti 9, 20133 - Milano
>> > tel.: +39-02-503-18304  fax: +39-02-503-18284
>> >
>> > SPIRITVS DVRISSIMA COQVIT
>> > _______________________________________________
>> > gmx-users mailing list
>> > gmx-users at gromacs.org
>> > http://www.gromacs.org/mailman/listinfo/gmx-users
>> > Please don't post (un)subscribe requests to the list. Use the
>> > www interface or send it to gmx-users-request at gromacs.org.
>>
>> _______________________________________________
>> gmx-users mailing list
>> gmx-users at gromacs.org
>> http://www.gromacs.org/mailman/listinfo/gmx-users
>> Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-request at gromacs.org.
>>
>--
>Groeten, David.
>________________________________________________________________________
>Dr. David van der Spoel, 	Dept. of Cell & Mol. Biology
>Husargatan 3, Box 596,  	75124 Uppsala, Sweden
>phone:	46 18 471 4205		fax: 46 18 511 755
>spoel at xray.bmc.uu.se	spoel at gromacs.org   http://xray.bmc.uu.se/~spoel
>++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
>_______________________________________________
>gmx-users mailing list
>gmx-users at gromacs.org
>http://www.gromacs.org/mailman/listinfo/gmx-users
>Please don't post (un)subscribe requests to the list. Use the
>www interface or send it to gmx-users-request at gromacs.org.
>

http://netwinsite.com/dbabble/
_______________________________________________
gmx-users mailing list
gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.