986

J. Chem. Theory Comput. 2005, 1, 986-993

Intrinsic Relative Stabilities of the Neutral Tautomers of Arginine Side-Chain Models Jan Norberg,*,† Nicolas Foloppe,‡ and Lennart Nilsson† Center for Structural Biochemistry, Department of Biosciences at NoVum, Karolinska Institutet, SE-141 57 Huddinge, Sweden Received December 18, 2004

Abstract: The specific protonation state of amino acids is crucial for the physicochemical properties of proteins and their biological functions. These protonation states influence, for instance, properties related to hydrogen bonding, solubility, and folding. pKa calculations for proteins are, therefore, important and require, in principle, a specification of the most stable protonated and deprotonated forms of each titratable group. This is complicated by the existence of multiple tautomers, like the five neutral tautomers of the guanidine moiety in arginine. In this study, the compounds N-methyl-guanidine and N-ethyl-guanidine were used to model the charged and all neutral protonation states of the arginine side chain. The relative stabilities of all five neutral tautomers were investigated systematically for the first time, using quantummechanical calculations. These relative stabilities were obtained in vacuo, water and chloroform, by combining the quantum-mechanical calculations with a continuum solvation model. The water model was used to represent arginines exposed to an aqueous solution, whereas the chloroform model has a polarity representative of a protein core or a membrane. This allowed determining the relative pKa’s associated with each neutral tautomer in these environments. A key result is that significant differences in stability are found between the neutral tautomers, in both water and chloroform. Some tautomers are consistently found to be the most stable. These findings will be helpful to refine pKa calculations in proteins.

Introduction The guanidine is found at the extremity of arginine side chains and is, therefore, an important biochemical building block. This guanidine functionality is usually assumed to be in the protonated guanidinium form because of its high intrinsic basicity, with an effective pKa of about 12.0 in aqueous solution at 298 K.1 There are several examples where arginine is found to play a key role in relation to the mechanism of action of a protein.2-5 For instance, the amino acid compositions of the interfaces of a large number of protein-protein6-8 and protein-nucleic acid9 complexes have * Corresponding author tel.: +46 8 608 9263, fax: +46 8 608 9290, e-mail: [email protected]. † Karolinska Institutet. ‡ Present address: Vernalis (R&D) Ltd., Granta Park, Abington, Cambridge CB1 6GB, U. K.

10.1021/ct049849m CCC: $30.25

been characterized, showing that arginine occurs frequently at the interface of the complexes.8,9 Other specific examples include the activation of G protein coupled receptors,10 the voltage gating of ion channels,11 a possible role for neutral arginine in porin channel selectivity,12 and the substrate binding of the human FLAP endonuclease-1 enzyme.13 Peptides containing arginine have an important role in proteomics projects, which use mass spectrometry for rapid and reliable identification of proteins.14-16 The fragmentation efficiency depends strongly on the protonation state of arginine. Despite being very basic, the guanidine functionality is included as a titratable site in methods that aim at determining the protonated state of proteins.17-20 These methods typically rely on the solution of the Poisson-Boltzmann equation in the framework of continuum electrostatics18,20 and calculate the pKa shift of a titratable group when it is transferred from © 2005 American Chemical Society

Published on Web 07/13/2005

Intrinsic Relative Stabilities

aqueous solution to the protein. The calculations involve both the protonated and the nonprotonated states of each titratable group in the protein environment. In these pKa calculations, one proton has to be removed from the guanidinium functionality of each arginine side chain. Details of the positioning of the polar hydrogens can have a significant impact on the pKa’s calculated in a protein environment,21-23 and there are, a priori, five acidic protons on the guanidinium group, corresponding to five different neutral guanidine tautomers. Because none of these five protons are rigorously chemically equivalent, it would be very useful to have a ranking of their relative intrinsic acidities. This information is, to our knowledge, not available, probably because it would be difficult to determine the corresponding microscopic pKa values experimentally. This, indeed, would require one to distinguish experimentally between each of these five neutral tautomers in aqueous solution. Alternatively, the guanidinium and each of the five neutral tautomers can be studied independently using a computational approach. Recent studies have shown that the relative pKa’s of a series of analogous compounds in aqueous solution can be calculated with reasonable accuracy, using high level ab initio quantum-chemistry calculations, in combination with a continuum model to represent the condensed phase dielectric properties.24-26 One study that tested several continuum models24 found a good agreement between experimental and calculated pKa’s when the aqueous environment was represented using the isodensity surface polarized continuum model (IPCM).27 The present report applies a similar protocol to rank the relative intrinsic acidities associated with the arginine side chain, modeled with two compounds, N-methylguanidine (Me-Gua) and N-ethyl-guanidine (Et-Gua), to test the effect of the length of the hydrocarbon side chain on the results. The relative acidities have been investigated in water, as well as in a medium with a low dielectric constant, to represent a protein core or a membrane environment. The results do indicate that there are significant stability differences between the neutral tautomers of an arginine side chain.

Computational Methodology The compounds used to model the arginine side chain were the protonated and neutral forms of Me-Gua (Figure 1) and Et-Gua (Figure 2). The protonated forms (guanidinium) of Me-Gua (Figure 1) and Et-Gua (Figure 2) are denoted MeGuaP and Et-GuaP, respectively. The five neutral tautomers are referred to as Me-GuaA (Et-GuaA) to Me-GuaE (EtGuaE) (Figure 1), with Me-GuaX/Y (Et-GuaX/Y) denoting any of the neutral forms. Although larger molecules, even arginine itself, would be tractable on modern computers, it is on purpose that the model compounds were kept as simple as possible. This is indeed needed to assess the true intrinsic properties of the guanidine group, without the complications of conformational issues that inevitably arise with larger compounds. Larger compounds also open the possibility of intramolecular interactions involving the guanidine moiety (e.g., folded structures), which would only obscure the results regarding the intrinsic preferences for the protonation states.

J. Chem. Theory Comput., Vol. 1, No. 5, 2005 987

Figure 1. Structures of the protonated (Me-GuaP) and neutral (Me-GuaA to Me-GuaE) forms of Me-Gua.

Figure 2. Gua. The series as tautomers

Structure of the protonated (Et-GuaP) form of Etsame naming scheme was used for the Et-Gua for the Me-Gua series; therefore, the neutral of Et-Gua are not shown explicitly.

Figure 3. Thermodynamics cycle used for deriving the relative pKa’s, using Me-Gua as an example.

The dissociation constant Ka of an acid HA into its conjugate base A- and a proton H+ is related to the Gibbs free energy change (∆G) of this dissociation reaction by ∆G ) -RT ln Ka ) G(H+) + G(A-) - G(HA)

(1)

where R is the molar gas constant and T is the temperature; our calculations used T ) 298 K. The absolute pKa for the dissociation reaction is then pKa )

1 ∆G 2.303RT

(2)

The absolute pKa is difficult to calculate, in particular, because of problems associated with the determination of G(H+).28-30 This study, therefore, focuses on the relative stabilities of the neutral tautomers of Me-Gua and Et-Gua. This can be investigated using a thermodynamic cycle (Figure 3), from which the pKa difference, ∆pKa, corre-

988 J. Chem. Theory Comput., Vol. 1, No. 5, 2005

sponding to the difference in acidity associated with two different neutral tautomers can be obtained. The cycle relates the acid dissociation in the condensed phase (∆Gcp) to its dissociation in the gas phase (∆Ggas) and to the solvation free energies (∆Gsolv) of all the chemical species involved. Accordingly, if the acid Me-GuaP dissociates to produce one of its conjugate bases Me-GuaX, the corresponding free energy difference and the associated microscopic absolute pKa in the condensed phase are given by ∆Gcp(Me-GuaX) ) 2.303RT pKa(Me-GuaX) ) ∆Gsolv(Me-GuaX) + ∆Gsolv(H+) ∆Gsolv(Me-GuaP) + ∆Ggas (Me-GuaX) (3) where ∆Gsolv(Me-GuaX) and ∆Gsolv(Me-GuaP) are the free energies of solvation of a neutral guanidine tautomer and the guanidinium, respectively. ∆Gsolv(H+) is the free energy of solvation for the proton, and ∆Gsolv(Me-GuaX) is the deprotonation free energy of the guanidinium in the gas phase. Subtracting eq 3, obtained with a neutral tautomer, Me-GuaY, from its counterpart for another neutral tautomer, Me-GuaX, yields ∆∆Gcp(Me-GuaX - Me-GuaY) ) ∆Gsolv(Me-GuaX) - ∆Gsolv(Me-GuaY) + ∆Ggas(Me-GuaX) - ∆Ggas(Me-GuaY) (4) A key point is that the free energies of solvation of the proton and the protonated guanidinium from eq 3 cancel when eq 4 is formed. Therefore, the relative pKa difference between two distinct neutral guanidine tautomers depends only on (i) the difference in free energy of solvation between these neutral tautomers [∆Gsolv(Me-GuaX) - ∆Gsolv(Me-GuaY)] and (ii) the difference in the gas-phase protonation free energies associated with these neutral tautomers [∆Ggas(MeGuaX) - ∆Ggas(Me-GuaY)]. The gas-phase protonation free energies were computed using ab initio quantum mechanics, and the solvation free energies have been calculated using a continuum solvation model.27 All calculations were performed using the Gaussian 98 suite of programs.31 The 6-31G(d,p) basis set was used for all compounds. Further, the 6-31+G(d,p) and 6-311+G(d,p) basis sets were also applied to test the influence of the basis set on selected calculations in vacuo. First, the geometries of all the compounds were optimized by energy minimization in the gas phase to the default tolerances of the program. The gas-phase energy-minimized geometries were then used directly for all the calculations in the condensed phase. This is expected to be a reasonable approximation, given that there is no torsional freedom within the model compounds that could be drastically affected by the presence of a solvent. Also, a recent study on similar compounds tested the effect of a solvent on the internal geometry, and this influence was found to be minor, as expected.32 The solvation free energies were obtained using the reaction-field IPCM,27 which previously led to satisfactory agreement with experimental results for the relative pKa’s in a series of pyridine derivatives.24 The solvation free energies were obtained at 298 K by calculating single-point

Norberg et al.

energies with the IPCM solvation model, for both water and chloroform. The dielectric constants of water and chloroform were set to 78.56 and 4.806, respectively.33 In the IPCM model, the isodensity cutoff value was 0.001 atomic units.24 All calculations were performed at three levels of theory: (i) Hartree-Fock (HF), (ii) second-order Møller-Plesset (MP2),34 and (iii) with the density functional theory using the B3LYP method.35 Gas-phase free energy differences were calculated by adding the vibrational, rotational, and translational corrections at 298 K to the differences in heat of formation obtained directly from Gaussian.31 This thermal free energy correction was determined from a normal-mode analysis of the previously geometry-optimized structures in vacuo at the HF and B3LYP levels of theory. The gas-phase proton affinity (PA) was computed as the change in the heat of formation and in the zero-point vibrational energy: PA ) Hf (Me-GuaX) - Hf (Me-GuaP) + ZPVE(Me-GuaX) - ZPVE(Me-GuaP) (5) where Hf(...) and ZPVE(...) are the heat of formation and the zero-point vibrational energy, respectively.

Results The properties of interest were determined for all compounds at the HF/6-31G(d,p), B3LYP/6-31G(d,p), and MP2/6-31G(d,p) levels in the gas phase and in water and at the MP2/ 6-31G(d,p) level in chloroform. For the Me-GuaD, Me-GuaE, and Et-GuaE model compounds, a tendency toward a more pyramidal conformation of the amino groups was observed, as has also been reported for the amino groups of guanidine and guanidinium cations.36 The HF and B3LYP results are presented for the sake of completeness, but our analysis focuses on the MP2 results. Overall, for a given model compound, the results for the different levels of theory were quite consistent and do reveal differences across the neutral tautomers. Basis Set Dependence and Bond Order Analysis. Given that the amino groups were found to be planar or pyramidal in different tautomers of methyl-guanidine with the 6-31G(d,p) basis set, the influence of larger basis sets on these geometries was tested. These calculations were done with the 6-31+G(d,p) and 6-311+G(d,p) basis sets and the B3LYP and MP2 levels of theory for Me-GuaB, Me-GuaD, and Me-GuaE. These were chosen because the Me-GuaB was planar at the MP2/6-31(d,p) level, whereas Me-GuaD and Me-GuaE were slightly pyramidal. Me-GuaB remained planar for both the larger basis sets. Also, Me-GuaD and Me-GuaE showed the same pyramidal tendency for the larger basis sets as they did for the 6-31G(d,p) basis set. Heats of formation computed for the three compounds using B3LYP and MP2 with the two larger basis sets were ranked in the same order as with the 6-31G(d,p) basis set. These differences in the heats of formation, for a given compound, computed with the three different basis sets were