|Home | About | Journals | Submit | Contact Us | Français|
The absolute pKa values of 24 representative amine compounds, including cocaine, nicotine, 10 neurotransmitters, and 12 anilines, in aqueous solution were calculated by performing first-principles electronic structure calculations that account for the solvent effects using four different solvation models, i.e. the surface and volume polarization for electrostatic interaction (SVPE) model, the standard polarizable continuum model (PCM), the integral equation formalism for the polarizable continuum model (IEFPCM), and the conductor-like screening solvation model (COSMO). Within the examined computational methods, the calculations using the SVPE model lead to the absolute pKa values with the smallest root-mean-square-deviation (RMSD) value (1.18). When the SVPE model was replaced by the PCM, IEFPCM, and COSMO, the RMSD value of the calculated absolute pKa values became 3.21, 2.72, and 3.08, respectively. All types of calculated pKa values linearly correlate with the experimental pKa values very well. With the empirical corrections using the linear correlation relationships, the theoretical pKa values are much closer to the corresponding experimental data and the RMSD values become 0.51 to 0.83. The smallest RMSD value (0.51) is also associated with the SVPE model. All of the results suggest that the first-principles electronic structure calculations using the SVPE model are a reliable approach to the pKa prediction for the amine compounds.
Amine compounds play a crucial role in numerous chemical and biological processes. Many biomolecules are amine compounds, e.g. neurotransmitters benzylamine (1), phenylethylamine (2), norephedrine (5), ephedrine (6), dopamine (7), norepinephrine (8), epinephrine (9), and serotonin (10) depicted in Figure 1. Cocaine (11) is also an amine compound and this widely abused drug can potently bind to various neuronal binding sites, including dopamine transporter (DAT), to block the neurotransmitter dopamine (7) reuptake. Cocaine is the most reinforcing of all drugs of abuse.1,2,3 The disastrous medical and social consequences of cocaine addiction have made the development of an anti-cocaine medication a high priority.4,5,6,7,8,9,10,11,12,13,14 The alkaloid nicotine (12), initially found in tobacco leaves, is the addictive amine compound that maintains tobacco smoking behavior.15,16,17 More than any other abused psychostimulant, nicotine addiction is the number one cause of preventable mortality and is responsible for over 4 million smoking-related deaths each year.18,19 Nicotine produces its effects on the central nervous system (CNS) by interacting with nicotinic acetylcholine receptors (nAChRs) that are essential for synaptic transmission. Neuronal nAChRs are members of a superfamily of ligand-gated ion channels, which modulate the function of many major neurotransmitter systems and thereby influence a broad range of brain functions, such as cognition, learning, and memory.19,20,21,22,23 A variety of agonists and antagonists of nAChRs have been reported in literature. Majority of the nAChR agonists and antagonists reported so far are also amine compounds.24,25 Many amine compounds that could be used as the nAChR agonists/antagonists are structural analogs of nicotine or cocaine. Selective nAChR agonists/antagonists have therapeutic potential in the treatment of Alzheimer’s disease, Parkinson’s disease, dyskinesias, Tourette’s syndrome, schizophrenia, attention deficit disorder, anxiety, and pain, as well as tobacco-use cessation agents.19,23,26,27,28,29,30,31,32,33
Each of these biologically important amine compounds has a protonable amine moiety and, therefore, has both the protonated and deprotonated molecular species in aqueous solution. It has been known that both the protonated and deprotonated molecular species could bind with a protein.34 Evaluation of the phenomenological binding of a protein with an amine compound should account for the microscopic binding of the protein with both the protonated and deprotonated molecular species. The phenomenological binding affinity is dependent on not only the microscopic binding affinities of the protonated and deprotonated molecular species, but also the pKa of the amine compound and the pH of the solution. Therefore, it is essential for determining the phenomenological binding affinity to reliably measure or predict pKa of the amine compounds.35
The pKa values of a wide variety of organic and inorganic compounds have been measured by different experimental methods,36,37 such as spectroscopy, potentiometry, conductimetry, and competitive reactions. Yet, the experimental pKa measurements of those classes of compounds always involve some artifice and approximations. A priori knowledge of pKa’s can be crucially important in situations such as those encountered in the development of biologically active molecules.37 Thus, it is of considerable interest to develop reliable computational protocols that can be used to calculate this property and complement the experimental techniques. In general, the computational modeling not only provides an alternative way38,39,40 to obtain pKa values beyond the experimental technical limitations, but also may help to better understand the relevant chemical/biological processes at molecular level.41–67 It is ambitious and difficult objective to accurately predict absolute pKa based on the first-principles electronic structure calculations. 68, 69 A particular challenge to the pKa prediction is the determination of the solvent effects in the first-principles electronic structure calculation. An additional challenge to the pKa prediction for the neurotransmitters and some other amine compounds is related to the different rotamers (see below) that coexist in solution.70
In order to identify a reliable first-principles approach to the pKa prediction for these types of biologically important amine compounds, we have performed various first-principles electronic structure calculations on both the protonated and deprotonated states of a series of representative amine compounds, including neurotransmitters, cocaine, nicotine, and anilines. In the first-principles electronic structure calculations, solvent effects on the Gibbs free energies and thus on the pKa values of the amine compounds in aqueous solution were accounted for by using four different solvation models. Comparison of various computational results with the corresponding experimental data allows us to identify the most promising computational approaches to the pKa prediction for the biologically important amine compounds.
For convenience, the protonated and deprotonated states of an amine compound are denoted by AH+ and A, respectively. We have the following dissociation process:
The equilibrium dissociation constant (Ka) for reaction (1) is determined by the corresponding change of the standard Gibbs free energy via
Depicted in Figure 2 is the thermodynamic cycle used to calculate the Gibbs free energy change (ΔGa) associated with reaction (1). According to this thermodynamic cycle, ΔGa is given by the following expression:
in which is the Gibbs free energy change corresponding to reaction (1) in the gas phase. ΔGsol(AH+), ΔGsol(A), and ΔGsol(H+) are the solvation free energies of the protonated state (AH+), the deprotonated state (A), and the proton (H+), respectively.
Specially, for each of the neurotransmitters (1 to 10), it has been known that about 95% of the molecular species in the equilibrium mixtures take mono-protonation state in aqueous solution. The first protonation occurs on the amine group of the side chain. For 2 to 10, each neurotransmitter mainly has three stable protonated conformations:71 one trans (TH+) and two gauche (G1H+ and G2H+) rotamers (see Figure 3). So we can presume that the concentration of AH+ (for 2 to 10) consists of all of the stable protonated rotamers in aqueous solution.
The different rotamers of a neurotransmitter (2, 3, … or 10) have different microscopic dissociation constants (denoted by KTH, KG1H, and KG2H for rotamers TH+, G1H+, and G2H+, respectively) and different concentration ratios in solution:
In Eqs.(8) to (11), ηTH, ηG1H, and ηG2H refer to the fractions of the concentration for rotamers TH+, G1H+, and G2H+, respectively. In order to calculate the total, phenomenological pKa in comparison with available experimental data, we need to know not only microscopic Ka values (i.e. KTH, KG1H, and KG2H) for all of the stable rotamers, but also their fractions of the concentration. ηTH, ηG1H, and ηG2H will be determined by their relative Gibbs free energies in solution. A combined use of Eqs.(4) to (6) gives
The phenomenological Ka and pKa can be expressed as follows:
In Eq.(15), KTH, KG1H, and KG2H are the microscopic Ka of rotamers TH+, G1H+, and G2H+, respectively.
Eq.(15) indicates that, to predict the phenomenological pKa of a neurotransmitter, we first need to determine the microscopic pKa for all of the rotamers and their concentration ratios. The microscopic pKa of all rotamers and their concentration ratios can be determined by using the calculated relative Gibbs free energies of all rotamers and their deprotonated structures in solution. The practical computational task becomes the calculations of Gibbs free energies of all relevant protonated and deprotonated structures of the amine compounds.
Geometries of all molecular species involved in this study were fully optimized by employing density functional theory (DFT) using Becke’s three-parameter hybrid exchange functional and the Lee-Yang-Parr correlation functional (B3LYP)72 with the 6-31+G* basis set.73 Harmonic vibrational frequencies were determined at the same B3LYP/6-31+G* level to confirm that the optimized structures were local minima on the potential energy surfaces and to evaluate the zero-point vibrational energy (ZPVE) and thermal corrections to the Gibbs free energy at T = 298.15 K. The geometries optimized at the B3LYP/6-31+G* level were used to carry out second-order Møller-Plesset (MP2) single-point energy calculations with the 6-311++G** basis set. All these electronic structure calculations in the gas phase were performed by using the Gaussian03 program.74
Four different self-consistent reaction field (SCRF)75 procedures were employed in the solvation calculations on the protonated and neutral states of amine compounds to evaluate their solvent shifts in aqueous solution. All of the solvation calculations were performed at the HF/6-31+G* level by using the geometries optimized at the B3LYP/6-31+G* level in the gas phase. The similar solvation protocol was used to successfully solve other interesting chemical and biological problems.76,77,78,79,80,81 The first SCRF procedure used in the present study is the recently developed surface and volume polarization for electrostatic interaction (SVPE) model 82 implemented in GAMESS program.83 The SVPE is also known as the fully polarizable continuum model (FPCM),76–81 because it fully accounts for both surface and volume polarization effects in the SCRF calculation. Since the solute cavity surface is defined as a solute electron charge isodensity contour determined self-consistently during the SVPE iteration process, the SVPE results (converged to the exact solution of the Poisson’s equation with a given numerical tolerance) depend only on the value of the contour for a given dielectric constant under a given quantum chemical calculation level.82a This single parameter value has been calibrated to be 0.001au.82b
The other three SCRF procedures employed in this study are the standard polarizable continuum model (PCM), 84 the integral equation formalism for the polarizable continuum model (IEFPCM), 85 and the conductor-like screening solvation model (COSMO)86 implemented in the Gaussian03 program. All PCM, IEFPCM, and COSMO calculations in this study were performed by using the default choices of the Gaussian03 program for the recommended standard parameters. The dielectric constant of solvent water used for all of the solvation calculations is 78.5.
The calculated total Gibbs free energy of a molecular species in aqueous solution was taken as the sum of the free energy calculated at the MP2/6-311++G**//B3LYP/6-31+G* (or at the B3LYP/6-31+G*//B3LYP/6-31+G* level) in the gas phase (including the thermal corrections to the Gibbs free energy at the B3LYP/6-31+G* level when T = 298.15 K and P = 1 atm) and the corresponding solvent shift determined by the SCRF calculation at the HF/6-31+G* level. Finally, prediction of the free energy change of reaction (1) also requires knowing the absolute solvation free energy of the proton (H+) in aqueous solution, ΔGsol298(H+), in addition to the free energies calculated for all of the molecular species mentioned above. Due to the inherent difficulty of measuring absolute solvation free energy of an ion, the reported “experimental” ΔGsol298(H+) values have a wide range from −252.6 to −264.1 kcal/mol.87 We recently calculated ΔGsol298(H+) by using a high-level, ab initio method of incorporating a hybrid supermolecule-continuum approach88,89,90,91 based on the same SVPE procedure used in the present study. ΔGsol298(H+) was predicted to be −262.4 kcal/mol88 which was used in the present study. The total Gibbs free energy of the proton is ΔGsol298(H+) plus the gas-phase free energy −6.3 kcal/mol (accounting for contributions from the proton translation to enthalpy and entropy) at 298.15 K.
The computers used to perform the calculations in this study are HP’s Superdome supercomputer (a shared-memory system with 256 processors) at University of Kentucky Center for Computational Sciences and IBM x335 34-processors Linux cluster and SGI multiprocessors Origin computers in our own lab.
Summarized in Tables 1 and and22 and depicted in Figures 4 to to77 are the pKa values determined by using the Gibbs free energies calculated with various methods, in comparison with the corresponding experimental data. The only difference between the results in Table 1 and those in Table 2 is the level of the gas phase energy calculations, i.e. MP2/6-311++G**//B3LYP/6-31+G* for Table 1 vs. B3LYP/6-31+G*//B3LYP/6-31+G* for Table 2. For the calculated pKa values in a same table, different columns refer to the different solvation models used to determine the solvent shifts of the Gibbs free energies (from the gas phase to the solution). Given in the tables are also values of the root-mean-square-deviation (RMSD) of the calculated absolute pKa values from the corresponding experimental data.
The RMSD values in the tables reveal that the pKa values calculated using the MP2/6-311++G**//B3LYP/6-31+G* energies are slightly better than those using the corresponding B3LYP/6-31+G*//B3LYP/6-31+G* energies. Based on the MP2/6-311++G**//B3LYP/6-31+G* energies, when the SVPE, PCM, IEFPCM, and COSMO models were used, the RMSD values are 1.18, 3.21, 2.72, and 3.08, respectively. The respective RMSD values become 1.30, 3.88, 3.30, and 3.67 when the MP2/6-311++G**//B3LYP/6-31+G* energies are replaced by the corresponding B3LYP/6-31+G*//B3LYP/6-31+G* energies. Further, no matter whether the MP2/6-311++G**//B3LYP/6-31+G* or B3LYP/6-31+G*//B3LYP/6-31+G* energies are used, the qualitative order of the RMSD values is the same, i.e. RMSD(SVPE) < RMSD(IEFPCM) < RMSD(COSMO) < RMSD(PCM). The smallest RMSD value is always associated with the SVPE model.
Although the RMSD values of the absolute pKa values calculated by using the PCM, IEFPCM, and COSMO models are as large as 2.72 to 3.88, all types of calculated pKa values linearly correlate with the experimental pKa values very well. When the MP2/6-311++G**//B3LYP/6-31+G* energies were used, we got the following linear correlation relationships:
The correlation coefficient (R) ranges from 0.9523 to 0.9824 for these linear correlation relationships. With the empirical corrections using Eqs.(16) to (23), the theoretical pKa values are much closer to the corresponding experimental data, as seen in Tables 1 and and2,2, and the RMSD values become 0.51 to 0.83. The largest R values (0.9824 or 0.9782) and the smallest RMSD values (0.51 or 0.56) are always associated with the SVPE model.
The absolute pKa values of cocaine, nicotine, 10 neurotransmitters, and 12 anilines in aqueous solution were calculated by performing first-principles electronic structure calculations that account for the solvent effects using four different solvation models. Within the examined computational methods, the first-principles electronic structure calculations using the SVPE model lead to the absolute pKa values with the smallest RMSD value (1.18). When the SVPE model was replaced by the PCM, IEFPCM, and COSMO, the RMSD value of the calculated absolute pKa values became 3.21, 2.72, and 3.08, respectively. All types of calculated pKa values linearly correlate with the experimental pKa values very well. The correlation coefficient (R) ranges from 0.9523 to 0.9824 for these linear correlation relationships. With the empirical corrections using the linear correlation relationships, the theoretical pKa values are much closer to the corresponding experimental data and the RMSD values become 0.51 to 0.83. The largest R value (0.9824) and the smallest RMSD value (0.51) are always associated with the SVPE model.
The research was supported in part by NIH (grant R01DA013930 to C.-G. Zhan), National Natural Science Foundation of China (No.20503008;No.20528201), and the Center for Computational Sciences at University of Kentucky.