By varying the availability of protons, i.e. the acidity of the media, the balance of the equilibrium can be shifted. This provides a measure of the ease of proton disassociation of a site in a compound, the disassociation (or ionization) constant pKa, defined by the equation:
Alternatively, the pKa of a site can be thought of the pH at which the protonated and deprotonated fractions are equal. If the pH is higher than the pKa, the site is mostly deprotonated, and if the pH is lower than the pKa, the site is mostly protonated.
To take an example, consider the drug aspirin (acetyl salicylic acid). The structure of aspirin as given in the Merck Index, WDI or Medchem is shown on the left, CC(=O)Oc1ccccc1C(=O)O. However, the pKa of the carboxylic acid group in aspirin is 3.49, which means that under nearly all physiological conditions, this group is almost entirely deprotonated. Hence, it is the physical properties of the molecule on the right, CC(=O)Oc1ccccc1C(=O)[O-], that is of most relevant to computational chemists.
|Blood (maternal umbilical)||7.25|
|Lacrimal fluid (tears)||7.4|
The fact that only the physiological ionization state is required greatly simplifies and increases the applicability of pKa prediction. There is little to be gained from knowing the pKa of an ionizable group if it is significantly greater than or less than pH 7.3. Similarly, a pKa that is close to 7.3, implies that a significant fraction of the molecule may be either protonated or deprotonated, and analysis should be applied to both forms.
Consider a docking calculation of the following molecule NSC577 from NCI95. It is unlikely that there are any accurate experimental values or Hammett equations of the effects of the aromatic ring systems on the ionizable iodine, but it is clear that the site is deprotonated (by about 16 log units!).
The currently implemented algorithm to determine physiological ionization is built upon a pKa predictor. The algorithm can start from either the input or neutral form of the molecule. The method then iteratively either protonates the most basic site or deprotonates the most acidic site until there are no remaining "acid" or "base" sites for the specified pH, or until a loop is encoutered. The pKa's of all sites must be recalculated on each iteration, so that pKa reflects the ionization state of the rest of the molecule. One pleasant feature of this approach is that it can determine whether a neutral molecule is predominantly zwitterionic.
The table below contains a number of atom types and associated values used in an atom type based predictor. The pKa is given for the both the conjugate acid and base. For convenience, the SMARTS-like constraint [!H0] is ommitted from all entries in the "Acid" column. The rows are grouped by element, and the first matching type is used. In the table below, no more than two shells of an atoms neighbour need be inspected to determine an approximate pKa.
|Conjugate Acid||Value||Conjugate Base||Model Compound|
The very generic nature of the atom types in the above table allow this approach to be applied rapidly to large datasets of diverse compounds.
This leads to the Hammett Equation for pKa:
Where Sigma is a constant assigned to a particular substituent, Rho is a constant assigned to the particular acid or base functional group and pKa0 is the pKa of the unsubstituted acid or base. For benzoic acids, pKa0 is 4.20 and Rho is defined to be 1.0.
In the original formulation, two constants need to assigned to each substituent, Sigmameta for meta-substitutions and Sigmapara for para-substitutions. This was soon extended to aliphatic systems by Taft, by introducing Sigmastar (also written Sigma*). Currently over 40 forms of sigma constant have been defined, but many of these corrolate extremely well with each other.
As a worked example, consider the pKa prediction of the compound shown below, 4-chloro-3,5-dimethylphenol.
The Hammett equation for phenol has pKa0 = 9.92 and Rho = 2.23. The Sigmameta for -CH3 is -0.06 and the Sigmapara is 0.24. Hence the predicted value of the pKa is 9.92 - 2.23*(0.24-0.06-0.06) = 9.70. This compares extremely well with the experimental value of 9.71.
Currently, a database of over 170 Hammett/Taft equations has been encoded as SMILES/SMARTS (including chirality) from the compilation given in Perrin and Serjeant's "silver book". Similarly, a database of sigma constants for several thousand functional groups has been compiled by Corwin Hansch and Al Leo, a early TDT version of which has been placed in Daylight "contrib" by John Bradshaw.
A major benefit of Hammet/Taft equations is their ability to handle special cases.
One common approach, is to extend the set of known substituents with sigma transmission equations. For example, Sigmastar of -CH2-R can be estimated as 0.41 * the Sigmastar for R. Similarly, -CH=CH- has transmission coefficient 0.51 and -C6H4- has coefficient 0.30. Similar schemes include the Exner-Fiedler method for aliphatic cycles, and the Dewar-Grisdale method for polyaromatic systems. However, this approach cannot be used when a terminal group in unparameterized.
A second approach is to use molecular orbital methods to estimate sigma values when there is no experimental data. Using the strong corrolation between sigmameta and sigmapara and charges calculated with MOPAC's AM1 hamiltonian, Peter Ertl of Novartis has been able to calculate sigma constants for over 80,000 organic functional groups. Fortunately, as described in my MUG97 talk, "Efficient Protein and Nucleic Acid Perception", there exist efficient algorithms for pattern matching large numbers of functional groups in constant time.
If the neutral state of the protein is taken as the reference state, its electrostatic energy in a given ionization state is:
where xi is 1 when the group i is ionized, and 0 when it is neutral; is +1 for bases, and 1 for acids; pKaintrinsic,i is the intrinsic pKa of group i; is the absolute value of the interaction energy of groups i and j; M is the number of ionizable groups in the protein; R is the gas constant; and T is absolute temperature.
The intrinsic pKa of group i is given by:
where is the difference between the free energy change of ionization for group i in the otherwise unionized protein and in a model compound of pKamodel,i. The intrinsic pKa thus represents the pKa of the group in the protein with all other titratable amino acids in the neutral state. Given values for and for all groups, the above equation provides the basis for the computation of pH-dependent properties of protein-ligand complexes.
Traditionally, this technique has only been used for determining the protonation states of amino acid sidechains in proetins, where the values of pKamodel,i are simply looked up in a small table. However, small molecule pKa prediction methods such as those described above could easily be incorporated, not only allowing the protein the affect the ionization state of the ligand but potentially the ligand affecting the ionization state of the protein.
|SMILES||Name||Exp. pKa||Calc. pKa|
The results show an impressive agreement using a very simple atom type based method. The R2 value is 0.800, and a standard deviation of 0.95. Far more accurate figures on this benchmark can be achieved using Hammett & Taft equations, where ACD Labs' ACD/pKa obtains an R2 value is 0.978, and a standard deviation of 0.39.
The following table shows the input and output structures of the physiological state algorithm above on the ten ligands used in Metaphorics' DockIt crunch.