Semi-continuum Solvation with SZMAP¶
Water Probe Orientations¶
SZMAP is a technology—based on semi-continuum Poisson-Boltzmann theory—for analyzing proteins and/or ligands using an explicit probe water molecule to estimate various thermodynamic properties.
The Poisson-Boltzmann (PB) continuum model of water [Grant-2001], [Gilson-1988], and [Sharp-1990] is simple, effective, and with OpenEye’s OEZap, fast. But, near strong charges or in regions of constrained geometry PB can be too good at solvating a surface, leading to over-estimates of solvation free energies [Sharp-1991]. PB also fails to model hydrogen bonds, “structured waters”, and other geometric considerations important to drug design.
One common approach to capturing these geometric aspects of solvation is to use all-atom simulations such as molecular dynamics with explicit solvent. For some molecular properties this may be a good way to go. But the very heavy sampling cost is often not repaid in improved results [Nicholls-2008], [Nicholls-2009].
For solvation, a more efficient approach is to probe the region of interest with a single explicit all-atom solvent molecule, embedded in and interacting with continuum solvent and with the molecular surface.
SZMAP is just such a semi-continuum solvation model. It uses OEZap to perform the PB calculations. Because the continuum solvent responds to the charges and geometry of the solvent probe water-water interactions are captured. Because the probe is constrained by its physical size and shape, structural details of solvent to solute interactions are captured.
Energy Calculations¶
By sampling a full range of orientations for the probe at any position, thermodynamic properties are calculated using statistical mechanics.
For the purposes of discussion, we’ll assume the solute being analyzed is a protein. The energy for each probe water orientation j is the sum of the ZAP protein:water Coulombic interaction energy, the protein and water desolvation penalties and the van der Waals energy. For speed, the protein desolvation term is not computed for each orientation, but is instead computed once per point using a spherical united-atom representation of the probe water.
The partition function Q is the usual sum:
The probability of each orientation is then:
The probe water orientation entropy difference from continuum water is:
Likewise, the probe water free energy difference from continuum water is:
And of course, the enthalpy difference is the sum:
Uncharged Water and Difference Energies¶
The values Szmap TK calculates for \(\Delta G\), \(\Delta H\) and \(T \Delta S\) represent differences between the explicit probe water and the continuum water. In other words, a positive \(\Delta G\) indicates that the continuum model over-estimates the cost of displacing water at that position. When trying to decide which modifications are most likely to increase binding affinity, we would prefer to compare the water probe to the ligand atom that displaces it, or at least a proxy for the ligand atom. In SZMAP, this proxy is the water probe without any of the energy terms that involve a charge on the water (leaving only \(P_{desolv} + VDW\)). Because these energy terms have already been determined, calculating the properties for this “uncharged” or “neutral” water requires no additional time, it’s free. The difference in free energy between standard and uncharged water is negative where standard water is more favored and positive where the uncharged water, our hydrophobic group proxy, is more favored.
This OEEnsemble_NeutralDiffDeltaG energy and
the OEEnsemble_VDW energy (which is subtracted out
of the difference calculation above) are the primary SZMAP terms
relevant to drug design.
Szmap TK can also calculate the difference between our probe and a vacuum bubble of the same size (defined as just the \(P_{desolv}\) term). Because protein desolvation is treated as constant over all probe orientations only the vacuum free energy and enthalpy differences are calculated, not the entropy differences.
Regions are said to clash if the sum of the van der Waals and interaction energies is greater than a cutoff (by default \(0 \, kcal \, mol^{-1}\)). Calculated ensemble energies in clashing regions are usually not very meaningful.
When interpreting these difference energies, it should be beneficial for a hydrophobic ligand atom to overlap a region where the neutral free energy difference is positive because this involves displacing water that is less favorable than the neutral probe. Overlapping a region where the neutral free energy difference is negative, on the other hand, requires the overlapping group to be at least as polar as water and making a similar level of interactions.
Stabilization¶
It is possible to use Szmap TK to calculate stabilization energies as a reaction energy
where bulk-solvent is defined as \(0 \, kcal \, mol^{-1}\), by subtracting neutral difference free energies from a point near the protein/ligand complex, the protein alone, and the ligand alone.
Negative stabilization free energy values indicate water is stabilized at that location in the complex, over and above the same location in the ligand and the apo pocket. Positive values are found where the water is destabilized in the complex. The stabilization or destabilization says nothing about the actual level of stability, so these energies need to be compared to others such as the apo or complex neutral difference free energy.