INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Myoglobin is one of the most widely studied proteins. Numerous x-ray structures with different myoglobin ligands and under different physical conditions were solved in the past. Here, we focus on a set of structures of carbonmonoxymyoglobin (MbCO), which are prepared at pH values of 4, 5, and 6 (Yang and Phillips, Jr., 1996). This set of structures shows that at low pH, His-64an important residue for O₂ specificity (Springer et al., 1994)swings out of the CO binding pocket. It is assumed that a protonation of His-64 goes along with this remarkable conformational change (Yang and Phillips, Jr., 1996; Johnson et al., 1996). In this work we study the protonation behavior of MbCO by the method of continuum electrostatics (Ullmann and Knapp, 1999), where first the electrostatic interactions are calculated by solving the Poisson-Boltzmann equation (Bashford and Karplus, 1990), and after that the large number of possible protonation states is sampled with a Monte Carlo method (Beroza et al., 1991). We have already successfully applied this method to investigate protonation and electron transfer processes in several other proteins (Rabenstein et al., 1998a,b, 2000; Vagedes et al., 2000).

The protonation behavior of myoglobin has been studied extensively by theoretical methods (Shire et al., 1975; Bashford et al., 1993; Yang and Honig, 1994; Sandberg and Edholm, 1999). However, in these studies conformational flexibility was not considered explicitly. To take into account the structural differences within the set of x-ray structures, the continuum electrostatics method must be extended to include an ensemble of multiple conformations (Ullmann and Knapp, 1999). There are numerous approaches to include multiple conformations in the calculations of protonation behavior (Buono et al., 1994; You and Bashford, 1995; Beroza and Case, 1996; Sham et al., 1997; Schaefer et al., 1997; Alexov and Gunner, 1997, 1999). However, these methods are not feasible here, because they account only for conformational changes of titratable groups, explicitly simulated water molecules, and/or hydrogens.

Here, a small number of protein structures is given with completely different sets of atomic coordinates. Hence, not only atoms of the backbone and non-titratable groups are displaced, but also the dielectric boundary is changed. In this work we present a new method to treat a given ensemble of a small number of arbitrarily different conformers of a protein. Please note that throughout this work the term "conformer" refers to a whole structure, not only to a single residue. This method is essentially a Monte Carlo version of the approach described by Rabenstein et al. (1998a) and Ullmann and Knapp (1999). Our aim is to show that our method generates consistent population probabilities of the different conformers and to calculate titration curves of the different titratable groups of MbCOespecially of His-64in agreement with experimental data. Our results will also have implications for the assignment of spectroscopically detected taxonomic substates (Johnson et al., 1996).

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Coordinates and parameters

In our calculations we used x-ray structures of sperm whale MbCO at pH values 4, 5, and 6 (Yang and Phillips, Jr., 1996). The PDB identifiers and resolutions of these structures are 1spe and 2.0 Å for pH 4, 1vxc and 1.7 Å for pH 5, and 1vxf and 1.7 Å for pH 6. For the residues 79-81, the structures 1spe and 1vxc contain two alternative conformations A and B, which we considered separately. Thus, we got five different structures called in the following 4a, 4b, 5a, 5b, and 6, with the digit corresponding to their respective pH value of preparation and the letter corresponding to the conformation of residues 79-81. We removed the sulfate ion and all water molecules from each structure. The influence of water was considered exclusively by a higher dielectric constant in cavities and outside of the protein because the orientation of the water molecules is not known, which makes their electrostatic effect uncertain (Rabenstein et al., 1998b). We considered as titratable groups the C-terminus, the N-terminus, the two propionates of the heme, and all aspartates, glutamates, arginines, lysines, tyrosines, and histidines, with the exception of His-93, which coordinates the heme iron and is therefore considered to be non-titrating. We used an all-atom model where all hydrogen atoms were treated explicitly with the exception of the hydrogen atom of protonated carboxyl groups. In this case the hydrogen atom was treated implicitly by distributing its charge symmetrically to the two oxygen atoms of the carboxyl group. The coordinates of hydrogen atoms were generated with CHARMM (Brooks et al., 1983). The positions of hydrogen atoms were energetically optimized using the CHARMM program with the CHARMM22 force field (MacKerell et al., 1992), while the heavy atom positions were fixed. For this optimization all titratable groups were in their standard protonation (i.e., aspartate, glutamate, heme-propionate, and the C-terminus unprotonated, arginine, histidine, lysine, tyrosine, and the N-terminus protonated). All atomic partial charges, including that of the heme, and the atomic radii were taken from the CHARMM22 force field (MacKerell et al., 1992). However, we needed to apply some modifications and endorsements: as mentioned above, the proton of protonated carboxyl groups is not represented explicitly, but simply by symmetrically distributing its charge to the atoms of the unprotonated carboxyl group (Table 1). We did this because we cannot know the exact binding site of the proton. We demonstrated in the past that this simplification works properly (Rabenstein et al., 1998a,b, 2000; Vagedes et al., 2000). A more detailed approach would be to introduce more than two states for a titratable group similar to the way it was done for histidines here (Bashford et al., 1993; Beroza and Case, 1996; Alexov and Gunner, 1997). For similar reasons, the deprotonation of Arg, Lys, and the N-terminus was represented by symmetrically removing a unit positive charge from the atoms of the titratable group (Table 1). For the charges of deprotonated tyrosine, which are not part of the CHARMM22 parameter set, we based our charges on quantum chemical calculations as described in Rabenstein et al. (1998b). Atomic partial charges that are not explicitly part of the CHARMM22 parameter set are listed in Table 1.

Titration of single conformers

First, we calculated the protonation pattern of all five structures separately. The theoretical background of the calculation of protonation patterns by solving the Poisson-Boltzmann equation (Bashford and Karplus, 1990) and sampling the possible protonation states with a Monte Carlo (MC) method (Beroza et al., 1991) is reviewed by Ullmann and Knapp (1999). We give here only an outline of the explanations, representing the most essential parts to understand the new methods applied in this work.

The energy Gⁿ of a protonation state n of a protein, which is characterized by the protonation state vector x_n = (x, x, ... , x), is given by Eq. 1 (Bashford and Karplus, 1990, 1991; Beroza et al., 1991; Yang et al., 1993; Beroza and Fredkin, 1996; Antosiewicz et al., 1996; Ullmann and Knapp, 1999).

(1)

Depending on the group µ, which can be protonated or unprotonated, x adopts the values 1 or 0, respectively. x₀ is the protonation state vector of the so-called reference state, where all titratable groups are in their uncharged protonation state, i.e., x is 1 if µ is an acid and 0 if µ is a base. The sums run over all N titratable groups. pK_intr,µ is the so-called intrinsic pK_a value of the group µ, and W_µ is the electrostatic interaction energy of the groups µ and  if both are charged. pK_intr and the symmetrical W-matrix are both calculated by solving the Poisson-Boltzmann equation on a grid (Warwicker and Watson, 1982). For doing this, we used the program MEAD (Bashford and Gerwert, 1992; Bashford, 1997). In our calculation, the protein is represented by a cavity with the dielectric constant _p = 4 containing the protein atoms with their atomic partial charges. In the remaining volume, the solvent is represented by setting the dielectric constant to _s = 78.5 and the ionic strength to I = 0.1 mol/l. For a discussion of the value for _p, see Warshel and Russel (1984), Warshel and Åqvist (1991), Honig and Nicholls (1995), Warshel et al. (1997), Rabenstein et al. (1998a), and Ullmann and Knapp (1999). For the solvation of the Poisson-Boltzmann equation, we performed grid focusing (Klapper et al., 1986) in two steps. Initially, we used an 80-Å cube with a 1.0 Å lattice spacing centered at the geometric center of the protein, followed by a 20-Å cube with a 0.25 Å lattice spacing centered at the considered titratable group. We used an ion exclusion layer of 2 Å and a solvent probe radius of 1.4 Å. The pK_a values used for the model compounds are listed in Table 2. We accounted for the  and  tautomers of histidine by considering histidines as double sites, as described by Bashford et al. (1993).

The protonation probability x_µ of the group µ can be calculated by evaluating the thermodynamic average over all 2^N possible protonation states n given by Eq. 2,

(2)

This thermodynamic average cannot be calculated exactly, since the number of possible protonation states of myoglobin (2^N with N = 73) is far too large. Instead, we sampled the set of protonation states using a Metropolis MC method (Beroza et al., 1991). This MC method is implemented in our program KARLSBERG (Rabenstein, 1999), which is freely available under the GNU public license from our webserver (http://lie.chemie.fu-berlin.de/karlsberg/). To increase sampling efficiency, pairs of titratable groups coupled by >3 pK_a units were treated by special double moves (Beroza et al., 1991), and triples of titratable groups, where at least two of the possible three pairs are coupled by >6 pK_a units, were treated by triple moves (Rabenstein et al., 1998a). The elementary MC step is the attempt to change the protonation state of a randomly chosen titratable group with subsequent evaluation of the Metropolis criterion. In the case of a double or triple move, the protonation state of all two or three titratable groups are changed before the Metropolis criterion is applied. An MC scan comprises as many MC steps as titratable single groups, doubles, and triples are available. We did one complete sampling for each 0.1 pH increment in the range from pH 3 to pH 7. For each sampling we performed 10,000 full scans, and after that 100,000 so-called reduced scans, where all titratable groups whose protonation probability during the full sampling differed from unity or zero by <10⁴ were fixed in their respective protonation state and excluded from further sampling. We calculated the statistical error of the protonation probability of a single group by evaluating its correlation function as described by Beroza et al. (1991). For each single group it was smaller than 0.003, and the statistical error of the total protonation of the whole myoglobin was always smaller than 0.05 protons. All calculations were done for T = 300 K.

Titration of a conformational ensemble

To account for conformational variability in computations of electrostatic energies, Eq. 1 must be modified, yielding Eq. 3, which is the energy G^n,l of the protonation state n in a certain conformation l (Ullmann and Knapp, 1999).

(3)

The interaction parameters pK and W depend now on the actual conformation l. The conformational reference energy G = G  G is the energy difference between an arbitrarily chosen, albeit fixed, reference conformation r and the actual conformation l. For the computation of G, the protein must be in its reference protonation state for both conformations r and l, i.e., all titratable groups are in their uncharged protonation state. If there are L conformations in the conformational ensemble, the total number of combined protonation and conformational states is L · 2^N (with N being the number of titratable groups). The thermodynamic average in Eq. 2 must now be evaluated for all L · 2^N states, yielding Eq. 4 (Rabenstein et al., 1998a).

(4)

The occupation probability _m of a certain conformation m within the conformational ensemble can also be calculated by a thermodynamic average:

(5)

Again, the thermodynamic average can only be calculated analytically for a small number of titratable groups and conformations as it was done by Rabenstein et al. (1998a). However, in this study the number of titratable groups is N = 73 and the number of conformations is L = 5, so that the number of states is 5 · 2⁷³. We also sampled this increased number of states by an MC method. We did this by introducing conformational MC moves in addition to the titration moves. However, in doing so, the conformational sampling proved to be not efficient enough. To overcome this problem, we applied parallel tempering (Geyer, 1991).

In parallel tempering the MC simulation is not performed for a single system, but for I noninteracting systems in parallel. Each system is a copy of the original single system with the only difference in temperature T. The systems are numbered from 1 to I in a way that T_i+1 > T_i (with 0  i < I). Titration and conformational moves are applied as before to each of the parallel systems independently. However, a third kind of move is introduced, a so-called tempering move. This move exchanges the protonation state x_i and the conformation l_i of a randomly selected system i(i < I) with the protonation state x_i and the conformation l_i+1 of system i + 1. A tempering move is accepted with the probability

(6)

where G_i is the energy G^n,l according to Eq. 3 for the current protonation state n and conformation l of system i. The parallel tempering procedure yields a canonical ensemble at the corresponding temperature for each parallel copy. Although the simulation of several ensembles is more costly, the tempering moves decrease the correlation within each ensemble, and hence decrease the statistical error.

In our MC simulations for the complete set of ensembles for all five myoglobin structures we applied parallel tempering at the temperatures 300 K, 400 K, 533 K, 711 K, and 948 K. For the evaluation of our results we considered only the simulation data at 300 K. We added one tempering move and 10 conformational moves per MC scan to the conventional titration moves. The simulations at higher temperatures are coupled to the 300 K simulation by the tempering moves such that the sampling efficiency of the 300 K simulation is increased as described before. For each sampling we performed 2000 full scans and 100,000 reduced scans, leading to a statistical error (at 300 K) of <0.02 for the protonation probability of single groups, of <0.1 protons for the total protonation of myoglobin, and of <0.02 for the occupancy of single conformations.

For the MC sampling of the conformational ensemble, reasonable values for G are required. These values can in principle be obtained from a molecular dynamics force field like CHARMM (Rabenstein et al., 1998a; Ullmann and Knapp, 1999). However, it is unlikely that this method will give reasonable results in our case. The structures used here are derived from crystallographic data and therefore inevitably involve small uncertainties in the atomic coordinates. Small structural changes within these uncertainties can dramatically change the energy value obtained from a force field. One could try to solve this problem by adjusting the structures according to the used force field, e.g., with a kind of energy minimization. It is, however, questionable in which way this should exactly be done without introducing an uncontrolled bias of the results. Therefore, we decided to use another method, where a modification of the given x-ray structures was not necessary.

First, we determined the G values only within the two alternative conformations in the x-ray structure from pH 4 (4a and 4b) and in the x-ray structure from pH 5 (5a and 5b), i.e., we simulated two-member ensembles consisting of 4a and 4b or of 5a and 5b. The occupancy of the two conformations was determined by the crystallographers to be 0.78 for 4a, 0.22 for 4b, 0.76 for 5a, and 0.24 for 5b. Our aim was to find G values that reproduce these occupancies, denoted as , in the multiconformational MC sampling. We achieved this by an iterative method. We first did MC simulations (one for 4a/4b, one for 5a/5b) where the initial value of G was set to zero. The pH values for these simulations were the same as those for the preparation of the x-ray structures, i.e., pH 4 for 4a/4b, pH 5 for 5a/5b. From the calculated occupancies of the conformations, a correction energy G was derived by Eq. 7, which was added to G.

(7)

If our initial MC sampling were not afflicted with a statistical error, the corrected values for G would be those we were looking for. However, there is a statistical error of , and the resulting error of G is especially large if is near to zero or unity, so we iteratively repeated the MC simulations with the corrected G values to get new values for G. We did this until the root-mean-square deviation (rmsd) of the calculated occupancies compared to the target occupancies was smaller than 0.01, which was fulfilled after one to three iterations, depending on the calculation. For convenience, we shifted the G values additively after each iteration so that the values for 4a and 5a were always zero. As G values, we obtained 6.78 kJ/mol for 4b and 0.11 kJ/mol for 5b.

In the next step, we determined G values for the sampling of all conformations together. Again we applied an iterative method as before. The problem is that values for the target occupancies are not directly available from experiment. Therefore, we assumed that, averaged over the pH range from 3 to 7, all three groups of structures, 4a/4b, 5a/5b, and 6, represent the complete ensemble of the infinite number of possible myoglobin conformations equally well. Hence, their occupancy integrated over this pH range should be equal (1/3). With this assumption, we could apply the iterative method as before with only a few modifications: the occupancy of a conformation is now calculated by integrating the calculated occupancy at a certain pH over the pH range from 3 to 7:

(8)

For this calculation, the occupancies of 4a and 4b (or 5a and 5b, respectively) were always added together. Their relative G values were fixed to the values given above. The initial guess for the G values was 100 kJ/mol for 4a and 5 kJ/mol for 5a. The value for structure 6 was fixed to be zero. The resulting G values were 98.29 kJ/mol for 4a, 91.51 kJ/mol for 4b, 2.05 kJ/mol for 5a, 2.16 kJ/mol for 5b, and 0.00 kJ/mol for 6.

Please note that the physical meaning of these energy values is rather limited. They result merely from the postulation that all three groups of structures, 4a/4b, 5a/5b, and 6, are populated equally averaged over the considered pH range, and that the population ratio of the structures 4a/4b or 5a/5, respectively, is given by the occupancies from the x-ray refinement data.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Population of the conformers

The calculated population probability of the conformers 4a, 4b, 5a, 5b, and 6 is plotted over the pH range from 3 to 7 in Fig. 2. Indeed, conformer 6 has its maximum population at the highest pH value and conformers 4a and 4b are most populated at the lowest pH value. Conformers 5a and 5b show a maximum of their population probability at intermediate pH values. For 5b, this maximum is located as expected at about pH 5, whereas for 5a it is shifted to a higher pH of ~6. This difference is remarkable because the two conformers differ only at the residues Lys-79, Gly-80, and His-81, and the rmsd is low (see Table 3). However, the shapes of the population curves of the corresponding conformer pair 4a and 4b differ less. Also, the differences between 4a/4b and the other conformers is not surprising because there is a larger rmsd of ~1.1 Å and the extreme displacement of the His-64 side chain. The conformational differences between 5a/5b and 6 are less pronounced, but nevertheless the shapes of their population curves are quite different.

View larger version (17K): [in this window] [in a new window]   FIGURE 1   Heme, CO ligand, and side chains of His-64 and His-93 from the three x-ray structures at pH 4, 5, and 6. Hydrogen atoms of the heme are not drawn for the sake of clarity. The main feature of the structural changes is the rotation of His-64 out of the CO binding pocket at low pH. According to our results, protonation of His-64 goes along with this rotation. The unprotonated His-64 prefers the  tautomeric form (drawn with the MOLSCRIPT program (Kraulis, 1991).

View larger version (19K): [in this window] [in a new window]   FIGURE 2   Calculated population probabilities of the different conformers in the pH range from 3 to 7. The number of each conformer corresponds to the pH value at which the respective x-ray structure was prepared.

View this table: [in this window] [in a new window]   TABLE 3   Root-mean-square deviations (rmsd) in angstroms for the five x-ray structures of MbCO considered in this work. Each pair of compared structures was aligned using the algorithm of Kabsch (1976)

In summary, we conclude that our method was able to calculate population probabilities of the MbCO conformers that are qualitatively correct and consistent with the experimental conditions of the x-ray structure determination. These population probabilities are sensitive to even small conformational changes.

Titration behavior of individual sites

Depending on the titratable site, the titration curves are more or less similar for the different conformations. For some sites, however, there are extreme differences. This is the case for His-64, which features the swing out of the CO binding pocket at low pH values. Fig. 3 shows the protonation probability of His-64 for different conformers. For the conformers 5a, 5b, and 6, His-64 is nearly unprotonated at all pH values. For the conformers 4a and 4b, where His-64 has swung out of the CO binding pocket, it is partially protonated, up to 80% at low pH values, but even at high pH values it is nearly 50%. If all conformers are considered together, the protonation of His-64 is closely correlated to the population probability of conformers 4a and 4b (Fig. 2) and thus with the rotation of His-64. This is in agreement with experimental findings (Yang and Phillips, Jr., 1996). The -proton of the protonated His-64 in the 4a/4b conformation is in hydrogen bond distance (1.8 Å) to the backbone oxygen of Asp-60. The distance to the carboxy group of Asp-60, which is according to our results negatively charged at all pH values, is 3.7 Å. The hydrogen bond and the electrostatic attraction between the positively charged His-64 and the negatively charged Asp-60 may be responsible for pulling the protonated His-64 out of the binding pocket. In the conformers 5a, 5b, and 6, His-64 is not involved in any strong hydrogen bond with the protein moiety or the heme. (However, hydrogen-bonding to water molecules in cavities of the protein is still possible. In the x-ray structure, water molecule 35 (structure 6) or water molecule 36 (structure 5a/b) is in hydrogen bond distance (3.1 Å or 2.7 Å, respectively) to N of His-64). In our calculations, we also determined the population probabilities for the two tautomers of histidine. Our result for His-64 is shown in Fig. 4. The unprotonated His-64 prefers the -tautomer, which is also reflected in Fig. 1.

View larger version (13K): [in this window] [in a new window]   FIGURE 3   Protonation probability of His-64 for the single conformers 4a and 6, and for the complete ensemble of conformers ("ensemble"). The protonation probabilities for the other conformers are not shown because they are very similar to those of conformer 4a (4b) or 6 (5a and 5b). His-64 is protonated to a significant amount only in conformers 4a and 4b. The protonation of His-64 for the whole conformational ensemble reflects the increasing population probability of the conformers 4a and 4b with decreasing pH (see also Fig. 2).

View larger version (14K): [in this window] [in a new window]   FIGURE 4   Population probability of the tautomers of His-64, considering the complete ensemble of conformers. The curve for the protonation probability is identical to the "ensemble" curve in Fig. 3.

The titratable residues Lys-79 and His-81 are directly involved in the conformational difference between 5a and 5b, or 4a and 4b, respectively. Their titration behavior may be a key to understanding the different shapes of the 5a and 5b population curves in Fig. 2. Lys-79 is completely protonated for all considered pH values and conformers. The titration behavior of His-81 is shown in Fig. 5. Although the titration curves for the conformers 4a and 4b are similar, the curves for 5a and 5b are much more different: the 5b curve resembles more the curves for 4a and 4b, whereas the 5a curve is similar to the curve for conformer 6. There are extreme differences between the a and b conformations in both cases (4a/4b and 5a/5b) for the neighboring residue His-82 (Fig. 6). The titration curve of 4a is similar to that of 5b, and the titration curve of 4b is similar to that of 5a, so the curves are swapped between the more populated a conformations and the less populated b conformations.

View larger version (18K): [in this window] [in a new window]   FIGURE 5   Protonation probability of His-81 for the different conformers and for the complete ensemble of conformers ("ensemble").

View larger version (17K): [in this window] [in a new window]   FIGURE 6   Protonation probability of His-82 for the different conformers and for the complete ensemble of conformers ("ensemble").

The strongly coupled pair His-24/His-119

As discussed by Bashford et al. (1993), the  nitrogen atoms of His-24 and His-119 are in hydrogen-bond distance and the two residues constitute a strongly coupled pair of titratable groups. During titration, they effectively share their -proton. Protonation of both residues at the same time nearly never occurs. In the unprotonated state, His-24 nearly completely adopts the -tautomeric form. The -tautomer of His-24 is nearly unpopulated. If and only if His-24 is protonated, His-119 adopts the  tautomeric form, so only three states are practically occurring: p, p, and (Fig. 7). A "p" means protonated, "" represents the -tautomer, and "" the -tautomer. The first letter denotes the protonation state of His-24, the second that of His-119. The states p and p are equivalent if the -proton is considered as shared. This behavior is illustrated in Fig. 8.

View larger version (29K): [in this window] [in a new window]   FIGURE 7   Structural representation of the strongly coupled pair His-24/His-119 with the three occurring protonation states p, p, and explained in text (drawn with the MOLSCRIPT program (Kraulis, 1991).

View larger version (16K): [in this window] [in a new window]   FIGURE 8   Population probabilities of the three possible states of the strongly coupled pair His-24/His-119 (see text). Because p and p are considered to be equivalent, their sum is also shown.

Comparison of calculated pK_a values with experimental values and other theoretical studies

There are experimentally determined pK_a values for most of the histidines and tyrosines in MbCO. The values for Tyr-103 and Tyr-151 were determined by Wilbur and Allerhand (1976) to be 10.3 and 10.5, respectively. These values are out of the pH range we considered here (pH 3-7). According to our calculations, both tyrosine residues are protonated over the whole considered pH range. This is in agreement with the experimental findings.

More interesting is the comparison of the histidine values. They are listed in Table 4. Our calculated values are compared to experimental values and calculated values from other models. In the null model, the pK_a values are simply set to the pK_a values of the model compounds in aqueous solution. Here, the  or  tautomer is chosen according to experimental results. In the uniform-80 model (Sandberg and Edholm, 1999), an electrostatic calculation is done with the dielectric constant uniformly set to  = 80 everywhere. SCPB denotes a single-conformer Poisson-Boltzmann calculation (Bashford et al., 1993). DaPDS (distance and position dependent screening) is a new method to calculate protonation states in proteins (Sandberg and Edholm, 1999). It is based on empirical screening functions (Warshel et al., 1984) and is much simpler and thus faster than methods where the Poison-Boltzmann equation must be solved. For each model we calculated the rmsd compared to the experimental values. In the calculation of the rmsd we included only histidine residues for which an experimental pK_a value in the range from 3 to 7 could be determined. In addition, we omitted the strongly coupled pair His-24 and His-119 due to their special titration behavior (s.a.), which does not allow assigning a well-defined pK_a value for each of these residues. That is also reflected in experimental findings (Bashford et al., 1993). According to our results, His-12 and His-81 have pK_a values >7, which is out of the range of our study. However, the point of half protonation is nearly reached at pH 7, so that we could extrapolate a pK_a value for these residues and include them in the rmsd calculation. The lowest rmsd value for the non-trivial models is reached with our method (Table 4). The low rmsd value of the uniform-80 model is not surprising because all the experimental pK_a values are not far from the model compound pK_a values, and the uniform-80 model trivially provides pK_a values close to the model compound values. Especially interesting are the pK_a values for His-64, which performs the swing out of the binding pocket. The SCPB calculation (Bashford et al., 1993), which used an x-ray structure prepared at intermediate pH (Kuriyan et al., 1986), where His-64 is situated inside the CO binding pocket, shows a very low pK_a value, as this is also implied by our calculation if only the conformers 6, 5a, or 5b are considered (Fig. 3). If we consider the ensemble of all five conformers, the calculated pK_a value of His-64 is in reasonable agreement with the experimental result. The DaPDS calculation (Sandberg and Edholm, 1999) yields a high pK_a value for His-64 of 6.30, although it uses the same x-ray structure as the SCPB calculation, where His-64 is situated inside the CO binding pocket and should be mostly unprotonated. This constitutes an example where simple methods like the uniform-80 model or empirical models like DaPDS yield, on the average, results in reasonable agreement with experimental findings, but do not provide an understanding of the underlying molecular mechanism. The calculation of pK_a values by solving the Poisson-Boltzmann equation may sometimes run into difficulties and yield deviating results, but is capable of unveiling what is going on in the more special situations, which may be crucial to understand protein function.

Bashford et al. (1993) reported a strong dependency of the Poisson-Boltzmann results on the used parameter set, i.e., atomic partial charges and atomic radii. They showed that choosing an appropriate parameter set can significantly improve the reliability of the results. The CHARMM22 parameter set, which was applied here, was not specifically designed for continuum electrostatics calculation, but rather for molecular dynamics at  = 1 and with an explicit solvent (if any). Possibly our results could be improved even further if the parameter set were tuned for continuum electrostatics calculations. However, our results presented here and also in another work (Vagedes et al., 2000) suggest that the CHARMM22 parameter set is already quite well suited for continuum electrostatics calculations.

Implications for the assignment of taxonomic substates

Three taxonomic substates of MbCO can be distinguished by infrared spectroscopy: A₀, A₁, and A₃. A₀ is more populated at low pH, whereas A₁ and A₃ are present at higher pH (Johnson et al., 1996). In consistency with our own results, the A₀ substate is usually assigned to the protonated His-64 rotated out of the CO binding pocket (for further discussion see Johnson et al. (1996), Müller et al. (1999), and references cited therein).

On the basis of this assignment, Müller et al. (1999) recently derived the pK_a values of His-64 and His-97 from their FTIR measurements of the taxonomic substates. These pK_a values are in good agreement with our calculated values (see Table 4).

In general, A₁ is the dominant substate at pH 5.7 in solution (Johnson et al., 1996). Ray et al. (1994) assigned A₁ to the -tautomer of His-64, and A₃ to the -tautomer. Jewsbury and Kitagawa (1994) criticized this assignment based on their MD simulations and suggested it is the other way around, i.e., A₁ is the -tautomer and A₃ is the -tautomer. Our results, however, support the model of Ray et al. (1994) because the -tautomer is preferred at higher pH (Fig. 4). However, it is also possible that the difference between A₁ and A₃ is not the protonation state of His-64, but a hydrogen bond to a solvent molecule. This is corroborated by the fact that the interconversion between A₁ and A₃ is very fast at room temperature, which cannot be explained by a simple tautomeric conversion. For more details see Müller et al. (1999).

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

We calculated the titration behavior of MbCO in the pH range from 3 to 7 by conventional electrostatic continuum methods with subsequent MC sampling of the protonation states. However, we considered not only one conformer, but an ensemble of five conformers derived from x-ray structures. These structures were solved for pH values 4, 5, and 6, so that the conformational ensemble used in our calculations represents the conformational changes of MbCO induced by pH changes. We extended the MC method by introducing conformational moves to sample the conformational ensemble in addition to the ensemble of protonation states. In our extended MC method the conformational ensemble may consist of a limited number of arbitrarily different conformations. The efficiency of the sampling was successfully increased by the application of parallel tempering (Geyer, 1991). As a result, we obtained population probabilities for the five conformers as a function of pH value. According to our results, the conformers corresponding to the x-ray data determined at pH 4 are most populated at low pH, whereas the conformer corresponding to pH 6 is most populated at high pH. The population probabilities of the conformers corresponding to pH 5 reach their maximum values at intermediate pH. Hence, our method was able to yield pH-dependent population probabilities of the conformers that are in qualitative agreement with experimental findings. The pK_a values calculated by sampling the conformational ensemble are in better agreement with experimental values than values calculated by a single-conformer Poisson-Boltzmann calculation (Bashford et al., 1993) and by the recently published DaPDS model (Sandberg and Edholm, 1999). We calculated the titration behavior of His-64, which features a remarkable rotation out of the CO binding pocket at low pH value, in agreement with experimental findings. Our results support the assignment of the taxonomic substates A₀ to the protonated rotated His-64, A₁ to the -tautomer, and A₃ to the -tautomer.