Yuri A. Kosinsky, Alexander S. Arseniev and Roman G. Efremov *
M.M. Shemyakin and Yu.A. Ovchinnikov Institute of Bioorganic Chemistry, Russian Academy of Sciences, Ul. Miklukho-Maklaya, 16/10, Moscow V-437, 117871 GSP, Russia.
* Corresponding author, phone: (7-095) 335 51 55; e-mail: efremov@nmr.ru.
The aim of this study is to employ electrostatic calculations to describe subtle differences in conformation and microenvironment of ionizable groups in neurotoxin II (NT) from Naja naja oxiana by fitting of their simulated and experimental pKa's values evaluated by measuring of chemical shifts of protons in the range 2-10 of pH. Spatial structure of NT which was solved earlier by 2D NMR (Golovanov et al., 1993) reveals unambigiuos conformations for charged side chains. Thus, caboxyl groups of Glu and Asp form transient hydrogen bonds with the backbone amide protons, and populations of these bonds were estimated from titration shift of respective amide protons. To select the conformers satisfying both NMR NOE distance restraints and observed pKa's values, the following computational protocol has been proposed: 1) conformational space of charged side chains was explored in detail via Monte Carlo simulation; 2) for each cluster of resulting structures, the values of pKa's were calculated using continuum electrostatics approximation; 3) further selection of the conformers was done by a comparison of experimental and calculated pKa's. The approach permits to reduce significantly the number of possible orientations of ionizable side chains obtained based on NOE data only. In addition, it might be employed to identification of charged groups in unfavorable environment - for which estimated pKa's deviate widely from those observed in the experiment. The method was tested on several ionizable groups of lysozyme, and the results reveal fairly good agreement with known X-ray data.
For many years, the calculation of pH-titration behavior of ionizable groups in proteins has provided a testing ground for theories of electrostatic interactions which are important for understanding of protein structure and function (reviewed by Warshel & Papazyan, 1998). Most of these calculations have been based on the assumption that the difference between the free energy of protonation of a particular group in protein and in small model compound can be treated in the framework of the quasi-macroscopic classical electrostatic model. The method of macroscopic electrostatics with atomic detail (MEAD) (Bashford & Karplus, 1990; Antosiewitcz et al., 1994) has been applied to calculate pH-titration behavior of ionizable groups in lysozyme and other proteins under assumption that there is little or no structural flexibility in or around the titrating groups. The conclusions reached through these calculations have demonstrated, however, that more accurate results are generally obtained by combining conformational flexibility data with macroscopic electrostatics approximation (You & Bashford, 1995).
Up to date, theoretical estimations of protein electrostatics are often made for "static" conformations - crystal structures or averaged NMR-derived models and, therefore, conformational lability of ionizable groups is not taken into consideration. On the other hand, knowledge of the details of conformation and charged state might have a paramount importance for delicate mechanisms of protein functioning (e.g., enzyme-substrate binding, catalysis, etc. (Warshel & Papazyan, 1998)). Although protein 3D structure and its mobility could be reconstructed based on NMR data, like NOE contacts and spin-spin coupling, in many cases this method fails to provide unambiguous assignment of orientations of ionizable groups. Solution to the problem could be achieved by simultaneous employment of structural constraints and pH-titration data derived by NMR spectroscopy, with theoretical methods.
In this paper we describe a computational approach to assessment of conformational mobility and electrostatic properties of ionizable groups in neurotoxin II (NT) from Naja naja oxiana. The method is illustrated by application to side chains of aspartic and glutamic acids, and includes the following steps: (i) exploring of conformational space of these residues in the three-dimensional (3D) structure of NT; (ii) theoretical estimation of their pKa's for the ensemble of NT structures that differ from one another in the regions of ionizable groups; (iii) refinement of conformations of ionizable groups by comparison of calculated pKa's with those evaluated by measuring of chemical shifts of protons in the pH-range 2-10. The approach was tested on a well-studied model - Glu residues in the X-ray structure of lysozyme. As we will demonstrate below, the method permits to reduce significantly the number of possible orientations of ionizable side chains obtained based on NOE data only.
Conformational sampling.
The spatial structure of NT which was solved earlier by 2D NMR NMR (Golovanov et al., 1993), reveals unambiguous conformations for charged side chains of Glu and Asp. To search for possible conformations of these side chains which comply with NOE inter-proton distance constraints, and to estimate their microenveronment, we employed Monte Carlo (MC) simulation in torsion angles space. This was done with the FANTOM (version 4.0) program (von Freyberg & Braun, 1991). During the simulations, only randomly selected dihedrals (except omega) of a given Asp or Glu and its two neighbours in the sequence were sampled, while the rest of the protein was kept fixed. The conformations accepted at each MC cycle were energy minimized via 150 steps of conjugate gradients algorithm. The ECEPP/2 potential function (Nemethy et al., 1983) and dielectric constant 
= 20 were used. Only non-identical low-energy conformers (differ in at least by 5o and 15o for backbone and side chain dihedrals, respectively) were retained for subsequent analysis. The same computational protocol was also applied to Glu7 and Glu35 of lysozyme. X-ray structure of lysozyme was taken in the Brrokhaven Protein Data Bank, PDB ( Shimanouchi, et al., 1977), entry 1LZT.
pKa's calculations.
pKa's were calculated using the MULTIMEAD (version 1.1.8) program (You & Bashford, 1995). Partial atomic charges and van der Waals radia were taken from the CHARMM (Brooks III, et al., 1983) library. The calculations were undertaken with the parameters imitating the conditions used for NMR studies of NT (Golovanov et al., 1993; Zagaitov et al., 1997): temperature 303 K, ionic strength 0.1 mM. Dielectric constants = 4 and = 80 were set for protein interior and solvent, respectively. Numerical solutions of the Poisson-Boltzman equation were estimated first on external grid (81 points, grid spacing 1.00 Å) with origin in the center of mass and all boundary points assigned potential 0.0, and next - on internal grid (81 points, grid spacing 0.25 Å) centered on the ionizable group.
pKa's of Glu in X-ray structure of lysozyme.
To test the method proposed in this study, we chose high-resolution X-ray structure of lysozyme because electrostatic properties of this protein were extensively investigated by a number of other research teams, and therefore, lysozyme might be considered as a very suitable reference model. Possible conformations of side chains of residues Glu7 and Glu35 were prepared using MC search as described above. Inspection of the conformers along with their total energies and calculated pKa's permits the following conclusions: (i) four clusters might be delineated in the 
1-2 space; (ii) for two the most populated clusters of Glu7, conformational energy and carboxylate pKa significantly vary - the calulated pKa is off by -0.2 ( 0.8 pH units from the measured one ( Table2). Energy variations inside these clusters could be explained by sensitivity of interaction of Glu7 with closely spaced side chain of Lys1 to changes of the dihedral angle 3.
Measured pKa of Glu35 is 6.1-6.3 (Bartik et al.,1994). This is considerably higher than the average pKa of carboxylates in proteins and could be induced by its heavily hydrophobic surrounding. In a case of Glu35, MC search revealed two clusters of conformers which differ by the values of 2 angle - 155o and -75o. For members of the first cluster (it corresponds also to the side chain conformation n the X-ray structure), estimated pKa is ~6.2, and thus, agrees well with the measured one ( Table2). On the other hand, calculated pKa for the second cluster is extremely high (~16). Such abnormal pKa is caused by proximity of COO- group of Glu35 to apolar side chain of Leu56 and its burial location with respect to the protein surface. Although the energies of conformers from the two clusters are quite similar, the large difference between calculated and experimental pKa's obtained for the second cluster makes it possible to conclude that its occurrence is highly improbable. Therefore, the results obtained for two residues Glu of lysozyme demonstrate that the constraints imposed by calculated values of pKa's permit discrimination of unproper conformations. Below we will consider application of the approach to refinement of NMR-derived spatial structure of NT.
Identification of transient hydrogen bonds in NT.
Analysis of pH-titration curves measured using NMR spectroscopy in the range 2(10 of pH enabled determination of pKa values of carboxylates of Glu and Asp residues, imidazole rings, and N-terminal amino group (Zagaitov et al., 1997). Using the method of Szyperski et al. (1994), prominent changes (corresponding to alterations of chemical shifts greater than 0.2 ppm) on the pH-titration curves for chemical shifts of the backbone amide protons of Glu2 (0.8ppm), Glu20 (0.3ppm), Asn22 (0.3ppm), Arg32 (0.4ppm), were assigned to occurence of a number of transient hydrogen bonds (H-bonds) between the amide protons and COO- groups of Asp and Glu. Assuming that the probability of H-bond is proportional to the amplitude of pH-dependent change of the chemical shift (Szyperski et al., 1994), we estimated approximate populations of the transient H-bonds in NT ( Table1).
Transient H-bond |
Population |
|---|---|
| Glu20 O - HN Glu20 |
20% |
| Glu20 O - HN Asn22 |
20% |
Asp30 O - HN Arg32 |
25% |
| Asp30 O - HN Arg32 |
50% |
Conformational and electrostatic properties of Glu and Asp residues in NT.
Glu2. The minimized conformations of Glu2 obtained in the result of MC simulation, were subdivided into six clusters, depending on values of their 1 and 2 angles. We found that the structures from the cluster-4 reveal significant deviation (~4 pH units) of calculated pKa from the experimental one. This can be induced by too close (and thus, unfavorable) contacts of carboxylate with the backbone atoms of Thr16 and Lys15. As a consequence, this cluster was removed from further analysis.
Glu20. Similar protocol employed for Glu20 leaded to delineation of four distinct clusters in the 1-2 space (Fig. 1). In spatial structures of clusters 1 and 2, COO- group of Glu20 is involved in transient H-bonds with HN:Glu20 and HN:Asn22, respectively (Fig. 3). As discussed above, existence of these H-bonds is supported by the experimental data (Table 1). Carboxylates of Glu20 in the structures attributed to clusters 3 and 4, do not participate in H-bonds, and differ only in their chi2 torsion angles. pKa values calculated for conformers from cluster-3 (coloured in magenta on (Fig. 3) are essentially (by ~1.5 pH units) higher than that measured in the experiment (Fig. 2). Most probably, this is caused by low accessibility of COO- group to solvent as well as by its hydrophobic environment created by residues Cys40 and Gly41. Taking into account these considerations, we concluded that the conformational space available for Glu20 is well represented by the clusters 1, 2, and 4. In addition, based on the experimental data, it is possible to estimate statistical weights of these clusters in the equilibrium mixture of the conformers of NT in solution. To the contrary, the structures from cluster-3 are rather less probable.
Asp30 and Asp57. Side chains of these residues are less flexible than those of Glu2 and Glu20. As follows from Table1, they could be involved in transient H-bonds whose appearance is related to change of 2 angles. Thus, COO- group of Asp57 characterized by 1 = -70o and 2 = -150o, forms H-bond with HN:Glu2, whereas that with 1 = -70o, 2 = 170o - does not. Fluctuations of calculated pKa for Asp30 probably, are determined by proximity of flexible charged side chains of His31 and Arg32. Side chain of Glu37 is buried inside the protein and, hence, demonstrates pKa which is systematically higher (by ~2.5 pH units) as compared to the experimentally measured value (Table 2). We suppose that this discrepancy could be explained by principal limitation of the continuum electrostatics models to calculate pKa's of buried residues.
Site |
pKa estimated |
pKa measured |
Difference |
|---|---|---|---|
| lysozyme (1LZT) | |||
| Glu7 (m) | 3.1+/-0.3 | 2.6-3.1 (2.85) | 0.23 |
| Glu35 (m) | 6.2+/-0.3 | 6.1-6.3 (6.2) | -0.03 |
| neurotoxin-II (NT) | |||
| Glu2 (m) | 4.1+/-0.7 | 3.7 | 0.42 |
| Glu20 (m)* | 4.0+/-0.6 | 3.7 | 0.34 |
| Glu20 (w)** | 4.1 | 0.40 | |
| Asp30 (m) | 3.1+/-0.3 | 3.2 | -0.06 |
| Asp30 (w) | 3.1 | 0.10 | |
| Glu37 (m) | 6.7+/-0.2 | 4.2 | 2.54 |
| Asp57 (m) | 3.6+/-0.3 | 2.7 | 0.88 |
| Asp57 (w) | 3.5 | 0.80 |
* (m) - mean estimated pKa values;**(w) - mean-weighted estimated pKa values (Table1).
Theoretical investigation of electrostatic interactions in NT enabled to reduce significantly the number of possible orientations of ionizable side chains of Glu2, Glu20 and Asp57 and demonstrated fairly good agreement with pKa's measured in the experiment. It is important to note that all studied conformers - even those rejected in the course of the analysis - satisfyed well the structural constraints derived from 2D NMR data. Therefore, the method permits refinement of the structures of ionizable groups in proteins as compared with standard NMR technique. The mean calculated pKa's obtained in this study for Asp and Glu residues of NT (after removing of "bad" clusters of conformers) along with the experimental pKa values, are summarized in Table2. It is seen that for seven ionizable groups, the root mean square deviation (rmsd) and correlation coeffitient between calculated and experimental pKa's are 1.02 and 0.79, respectively. (For comparison, "null" model gives the values 1.10 and 0.48, respectively.) These data agree well with the results obtained for high-resolution X-ray structures of Asp and Glu resuides in lysozyme.
We thank Drs. W. Braun and D. Bashford for providing us with the FANTOM and MULTIMEAD programs, respectively. This work was supported in part by the Russian Foundation for Basic Research (grant 98-04-48823).
Fig. 1. Clustering of conformers of Glu20 of NT obtained in the result of Monte Carlo simulation. Each point on the (1-(2 map represents an energy minimized structure.
Fig. 2. Difference between calculated and measured pKa's for the clusters of conformers of Glu20 of NT.
Fig. 3. Skeletal drawing of representative conformers from different clusters of Glu20. Structures from clusters 1-4 are coloured in orange, green, magenta, blue, respectively.
