Solvation energy of a large atom cluster: continuous solvation energy test - II

As it has been already stated here, binding energy calculation of a small molecule to a large protein poses a difficult problem: a method for molecular electrostatic energy calculation should work well both for the protein ligand complex, the protein and the ligand at infinite separation. The protein and the complex are large molecules, whereas the ligand is, by definition, small.

Not every computational approach for the solvation energy calculation is fit for the job though. To elucidate the nature of the problems at hand we performed the following model calculation:
- we prepared a spherical "protein" of a large (but realistic) radius
- we placed a single-atom ligand with a charge at a given distance from the "protein" center (see the Figure)
- we calculated the solvation energy of the system as a function of the ligand distance both when the protein is neutral and charged (in the latter case the protein charge is opposite to that of the "ligand")

We used four different methods for the electrostatic contribution to the solvation energy calculation: Poisson equation solver (in its surface electrostatic incarnation, blue) QUANTUM MGB (cyan) and the two "classic" GB methods, based on the Coulomb approximation: HCT (magenta) and AGBNP (yellow).

The first Figure, corresponding to an overall electrically neutral cluster, shows absolute deficiency of HCT approach deep enough inside the "protein". The problem is caused by unrealistic assumptions with regard to the overlap integrals calculations is occurs pretty frequently in realistic proteins. AGBNP method represents one of the latest GB approaches and is practically free of these difficulties. However, AGBNP is based on Coulomb approximation and thus fails to recover correct behavior of the solvation energy close to the "protein" boundary: AGBNP energy is off by a large number from both QMGB and the exact solution. QMGB and Poisson solutions agree very well everywhere!

The last Figure shows the same calculation for a charged model "protein-ligand" complex. Once again, HCT fails entirely, AGBNP does not work properly at the "protein" boundary and both Poisson solver and QMGB agree very well, though QMGB is about one order of magnitude faster than the Poisson solver!

Practically all this means that QMGB represents a fast and accurate approximation to the Poisson equation solution. QMGB approach does not rely on Coulomb approximation and is shown to work both for small molecules and large molecular clusters involving molecules of very different sizes. Therefore, with QMGB one can find a single transferable set of GB parameters capable of describing large and small molecules on the same footing.

Solvation energy of a diatomic molecule: continuous solvation energy test - I

Diatomic molecule is the simplest example of a realistic solvation energy calculation. Indeed, any reasonable solvation energy model gives exact value for a single atom.

Depending on the radii of the atoms involved the solvation energy of a pair may be a very good test of a solvation energy model and transferability of its parameters. In what follows we show the results of our Modified GB (MGB) approach for the test system. The graph below represents the solvation energies for similar and opposite charges pairs.

The upper (blue) solid curve represents the atom pair with opposite charges, whereas the lower (red) curve corresponds to the atoms with similar charges. First of all, the behavior of the two energies is easy to understand. At infinite separation both curves saturate at -0.125 (which is the Born solvation energy of the two charges of 0.5 and bare radii 2.). If the total charge is 0 (the opposite charges case, blue curve), at r=0 we have Es=0 as well. If the total charge is 2x0.5=1 (the red curve), then at r=0 we have Es=-0.25, as it should be for a combined charge within the sphere of radius 2.

Although the asymptotic values are OK, this does not mean the whole curve is fine. To compare our approach with true electrostatic we performed the calculation of the model system solving the Poisson equation as well as by two "classic" GB models (that of HCT and AGNP). The results for a diatomic molecule with zero total charge are represented on the last graph.

The electrostatic part of the solvation energy corresponds to the blue curve of the previous graph and is calculated either by a (surface-electrostatic) Poisson equation solver (blue), QUANTUM's MGB (cyan), AGBNP (yellow) and HCT GB model (yellow). As it is clear from here, all the approaches give very similar results for the "small" molecule and are practically indistinguishable. Indeed, it is well known that practically any sort of GB approximation gives good results for solvation energies of small molecules.

The difference between QMGB method and "classic" GB approaches and its relation to the exact solution becomes more obvious if we consider a charged diatomic molecule (a molecular ion with total charge, say, 1 placed on one of the atoms). The exact (blue) and QMGB (cyan), once again, are both in agreement with each other, whereas both "classic" GB approaches, HCT and AGBNP fail to recover correct asymptotic value at zero inter-atomic separation. The latter difference between GB solutions and the exact value of the solvation energy is not important for small molecules (low atom density) but is extremely important for ligand binding calculations (to be explained).

Self-consistent solvation energy contribution calculation for protein-ligand complexes

Solvation energy is a major contribution to a ligand binding energy and is the interaction pretty much responsible for binding selectivity. Actual calculation of the solvation energy requires a method valid both for small molecule ligands and large proteins (and protein-ligand complexes).


Calculation of the electrostatic contributions for the binding energies in a continuous solvation energy approach imposes different problems for large and small molecules. Normally people use some kind of Generalized Born (GB) approximation. The latter is only exact for a charge in the center of a spherical cavity and thus can only be valid for a small molecule with most of the charge located within a few atoms.

If the molecule of interest is large, most of the charges are close to the molecular surface, instead. GB approximation in its most commonly accepted form fails next to a molecular surface: the Born radius is missed by a factor of 2. This means that there can be no "classic" GB model working good both for small and large molecules!

Binding affinity calculation requires calculation of differences between the energy of the complex (a large molecule) and the energies of the protein (another large molecule) and the ligand (a small molecule) at infinite separation.

If a GB model is made working by careful adjusting of "bare" Born radii to fit experimental IC50 of complexes, a good sanity check would require reproduction of experimentally known solvation energies of small molecules and ions (and the other way around). The two graphs in this post show, that this is indeed possible. A relatively large error in the small molecules solvation energies shows that although the resulting model is reasonable, the obtained GB parameters are only quantitatively transferable between large and small molecules.