A quantum-chemical picture of hemoglobin affinity
- R. E. Alcantara†,‡,
- C. Xu§,
- T. G. Spiro§, and
- V. Guallar‡,¶
- †Center for Computational Biology, Washington University, 700 South Euclid, St. Louis, MO 63110;
- ‡Life Sciences, Barcelona Supercomputing Center, c/ Jordi Girona 31, 08034 Barcelona, Spain; and
- §Department of Chemistry, Princeton University, Princeton, NJ 08544
-
Edited by Harry B. Gray, California Institute of Technology, Pasadena, CA, and approved September 25, 2007 (received for review June 28, 2007)
Abstract
Understanding the molecular mechanism of hemoglobin cooperativity remains an enduring challenge. Protein forces that control ligand affinity are not directly accessible by experiment. We demonstrate that computational quantum mechanics/molecular mechanics methods can provide reasonable values of ligand binding energies in Hb, and of their dependence on allostery. About 40% of the binding energy differences between the relaxed state and tense state quaternary structures result from strain induced in the heme and its ligands, especially in one of the pyrrole rings. The proximal histidine also contributes significantly, in particular, in the α-chains. The remaining energy difference resides in protein contacts, involving residues responsible for locking the quaternary changes. In the α-chains, the most important contacts involve the FG corner, at the “hinge” region of the α1β2 quaternary interface. The energy differences are spread more evenly among the β-chain residues, suggesting greater flexibility for the β- than for the α-chains along the quaternary transition. Despite this chain differentiation, the chains contribute equally to the relaxed substitute state energy difference. Thus, nature has evolved a symmetric response to the quaternary structure change, which is a requirement for maximum cooperativity, via different mechanisms for the two kinds of chains.
Hemoglobin has been a paradigm for our understanding of multisubunit proteins, beginning with the pioneering work of Adair (1) and then with the elucidation of the oxygenated and deoxygenated crystal structures of Hb by Perutz (2), some four decades ago. Interest in Hb continues to the present, with new experimental and theoretical methods aimed at characterizing intermediate ligation states of the molecule (3–5) and at gaining a deeper understanding of its allosteric transition (6–11).
The Hb molecule is composed of four chains that contain a heme group capable of reversibly binding CO or O2. The molecule is arranged as a dimer of dimers, where each dimer is made up of an α- and a β-chain (denoted as α1β1 or α2β2) with a hydrophobic intradimer interface. Adair observed that Hb binds oxygen in a cooperative manner, successive ligation events occurring with higher affinity. For some time it was believed that such cooperativity could be explained by direct interaction of the closely packed heme groups. However the 1960s saw the appearance of the first allosteric models for Hb cooperativity. In particular, the Monod–Wyman–Changeux (12) model postulated the existence of two states of the molecule: the relaxed (R) state with high ligand affinity and the tense (T) state with low affinity. Ligand binding preferentially stabilized the R form and led to cooperativity. Early Hb crystallization attempts revealed that deoxy crystals would shatter when exposed to oxygen, which suggested that a quaternary rearrangement was associated with the T-R transition of the molecule.
Crystal structures for oxy- and deoxy-Hb provided an atomic description of the structural changes that accompany oxygenation (2, 13). These structures were identified with the two functional states, R (oxy) and T (deoxy), of the Monod–Wyman–Changeux model. The most obvious difference involved the quaternary structures; one of the αβ-dimers is rotated ≈15° relative to the other. The interdimer interfaces are rearranged, whereas the intradimer interface remains largely unchanged. Two critical contact regions in the α1β2 interface were noted: the “hinge,” where rotation about the contacts between the C-helix of the β-chain and the FG corner of the α-chain reorient the residue side chains; and the “switch,” where the contacts between the C-helix of the α-chain and the FG corner of the β-chain are translated to an altered interdigitation register (2). Perutz also noted several salt bridges in the T structure that are broken in the R structure (13).
In addition, ligation-induced displacements of the heme iron and the proximal histidine ligand were detected in the crystal structures. These changes were in line with chemical expectations from porphyrin model systems, that pentacoordination of the heme in deoxyHb favored the high-spin state (quintet), with the Fe(II) displaced from the heme plane, whereas addition of a sixth ligand induced movement of the Fe into the plane and favored the low-spin state (singlet). Movement of the iron and proximal histidine are at the core of the Perutz model, constituting both the allosteric trigger and actuator for oxygen ligation. All of these elements taken together, reorientation of the dimers, breaking of the salt bridges, and proximal histidine displacement, constitute the classical Perutz model of Hb allostery (2). The Perutz model was supported by subsequent experimental and computational work (14). Some of the first Hb computational studies were carried out by Gelin and Karplus (15) and by Warshel (16) and later by Warshel and Weiss (17), who extended the original model. Karplus coined the term “allosteric core” to refer to those elements of the protein that carry stress from the ligation event (principally residues of the F-helix). Although great progress has been made in understanding the structural basis of affinity modulation in Hb, many details are not well understood.
This study examines the influence of protein structure on the energetics of the quintet and singlet states of the heme, as a ligand is bound or dissociated (Fig. 1), and how this influence is expressed in the transition from the T to the R state. Quantum-chemical methods have already been successfully used to study members of the globin family (18–20). In this article we employ a quantum mechanics/molecular mechanics (QM/MM) method to look at this coupling of conformation and spin change. QM/MM methods allow a quantum mechanical treatment for a subsection of a system, whereas the rest is treated with a classical description. This partitioning allows a high level of detail where it is required, for example, the ligand binding site, and a low-level description for the rest of the system. In the case of Hb, the quantum description is localized at the heme, whereas the rest of the protein is treated classically.
Potential energy profiles as a function of iron–CO distance. The two curves denote the two electronic states that are favored when the iron is hexacoordinated (singlet), or when it is pentacoordinated (quintet). The modulation of ΔE by the protein allows for changes in affinity of the system.
The results from our study identify those structural elements that control affinity. They offer a basis for understanding the effects of ligand binding on key structural elements that contribute to the allosteric transition (which should be explored in future work).
Results and Discussion
Affinity of the R and T States.
Table 1 lists computed energy differences between quintet and singlet states (ΔE) when CO (a much-studied surrogate for O2) is dissociated from HbCO (ΔE R) or bound to deoxyHb (ΔE T). These were calculated by in silico treatment (see Methods) of the highest-resolution crystal structures available for the two forms of Hb (21). To dissociate HbCO, we stretched the Fe–C bond until the CO dissociated and adopted a position adjacent to and parallel to the heme, the geometry observed in photodissociation experiments of myoglobin (Mb) (22) and Hb (23, 24). In deoxyHb a CO molecule was allowed to diffuse into the heme pocket and then bind to the heme. In the subunit subjected to these ligation changes, the heme and protein were allowed to relax via successive rounds of QM/MM optimization, whereas the remaining three subunits were frozen, to maintain the quaternary structure, R or T. The values presented in Table 1 are averages from five independent simulations, as indicated in Methods. With this procedure we obtain local minima describing affinity snapshots of the fully ligated and unligated crystal structures. Further conformational sampling is possible, for example, through extensive side-chain orientational sampling, as shown in our recent study on Mb (20). Such conformational sampling, however, could introduce mixing of the R and T structures and complicate the connection to heme electronic energies (because of large variations in the classical component of the QM/MM energy).
Total and partitioned ligation energy (kcal/mol)
The ΔE values are all positive (Table 1), indicating greater stability for the ligated singlet state, as expected. Moreover, the magnitudes, 9–19 kcal/mol, are in reasonable accord with the calorimetrically measured heats released on CO binding, 22 and 24 kcal/mol for isolated α- and β-chains (25). The measured values include transfer terms from solvent to the protein interior, whereas the computed ΔE values do not. Correction for these effects would improve the agreement with computed values.
The computations indicate higher stabilities for ligated α- than β-chains within the tetramer, by ≈2–3 kcal/mol in both the R and T structures. Experimentally, it is difficult to determine binding affinities for the two chains separately within the tetramer, and attempts to do so have yielded a range of results. The most thorough study appears to be that of Unzai et al. (26), who measured rates of ligation and dissociation for several hybrid Hbs, constructed to maintain the R or T quaternary structure. They concluded that for both CO and O2, the α- and β-chain affinities were similar, within a factor of 2 (≈0.4 kcal/mol of free energy), in both the R and T states. However, free-energy differences could entail significant entropic contributions, which are not included in the computation, in addition to the transfer terms mentioned above.
Despite the higher computed ΔE for α- than for β-chains, the R-T differences (ΔΔE; Table 1) are remarkably similar. In this respect, the computation is in good accord with the conclusions of Unzai et al. (26). The magnitudes of ΔΔE, ≈7 kcal/mol, exceed the experimental free-energy difference, ≈4.5 kcal/mol, but again entropy differences come into the comparison. The role of proton release (Bohr effect) also is important in the R-T energetics (25), and these are not considered in the computation.
Spurious effects caused by the constraining scheme, that is, freezing all chains but the ligating one, were evaluated by repeating some of the calculations without constraints. Results from these simulations showed that the constraints have a negligible effect; affinity values for these simulations were the same (within measurement error) of those previously found.
We conclude that the computed ΔE values are in the right ballpark and that they correctly point to similar affinity differences between the R and T structures for both chains.
The in silico procedure allows us to evaluate separate contributions to the computed ΔE values from the heme ligation itself and from consequent heme–protein and protein–protein interactions (Table 1). Heme ligation is the dominant contribution, as expected, but the other interaction terms are not insignificant. When R-T differences are taken, we find that the other interaction terms contribute ≈60% to the ΔΔE values; the remainder are due to effects on the heme itself.
Protein Contribution to Affinity Change.
To gain a deeper structural understanding, we looked further into the factors controlling ΔE. Initially, we examined the effect of electrostatics by transferring the heme from a high dielectric medium (solvent) into the electrostatic field of the protein. The effect was to stabilize the quintet state, by ≈2 kcal/mol, suggesting that electrostatics plays a role in the generally observed affinity reduction in heme proteins relative to heme in solution (27). This effect was essentially the same for both chains, in both the R and T structures. However, the early computational study by Warshel and Weiss (17) did find a differential electrostatic contribution of ≈1 kcal/mol when O2 was the ligand. Bound O2 is well known to experience a larger electrostatic attraction to the protein than bound CO does. Accounting for protein discrimination between the two ligands (28, 29), it is plausible that a small α-β difference in electrostatic energy would appear with the substitution of O2 for CO.
Next, we examined the protein–protein interactions by tabulating the contribution of each residue to the computed ΔΔE. Fig. 2 shows the disposition within the T structure of the α-chain residues having the largest positive contributions (values listed in Fig. 2 Inset; there are numerous small positive and negative contributions), while Fig. 3 does the same for the β-chain. Among the top seven contributors (>0.5 kcal/mol) in the α-chain, four (R92, K90, H89, and L86) are at the FG corner, which has long been identified as part of the allosteric core (15). The first two form intersubunit H bonds (Arg-92α-Glu-43β and Leu-91α-Arg-40β) as part of the hinge quaternary contact and are broken in the R state. Together they contribute 2.9 kcal/mol to ΔΔE, consistent with mutational evidence that the hinge is critical to the stability of the T state (30). On the other hand, His-89 and Leu-86 form an H-bond chain that stabilizes the proximal imidazole ligand (His-87), which is also an important contributor to the R-T energetics (see below). The distal histidine (His-58) is required to make a significant rearrangement to accommodate the ligand in the T state. It is buttressed in the T state by the adjacent Lys-61α, another major contributor to ΔΔE, which additionally forms a salt bridge with a heme propionate. Finally, Leu-136α is near the C terminus (Arg-145α), which is more mobile in the R than the T structure (31), where it is restrained by a quaternary salt bridge. Removal of residues near the C terminus has long been known to destabilize the T state (32, 33).
Residues with the largest contribution to the decrease in affinity ΔΔE of the α-chain shown with a stick representation. The heme, proximal histidine, and indirectly involved residues appear in a darker shade. (Inset) The value of the residue contribution (kcal/mol) is shown. Important interactions are shown with a dashed line.
Residues with the largest contribution to the decrease in affinity ΔΔE of the β-chain shown in a stick representation. Heme, proximal histidine, and indirectly involved residues appear with a darker shade. (Inset) Shown is the value of the residue contribution (kcal/mol). An important interaction is represented by a dashed line.
Turning to the β-chain (Fig. 3) we find that the ΔΔE contributions are more widely distributed among the residues than in the α-chain; none of the β-chain residues contribute >0.5 kcal/mol. Two of the top positive contributors, Glu-90β and Leu-88β, are on the F-helix, and support H-bond rearrangements (involving the proximal His-92 ligand in the case of Leu-88β) that increase mobility near the C terminus, His-146β, which also contributes significantly to ΔΔE. As in the α-chain, the C terminus is more mobile in the R than in the T structure, where another quaternary salt bridge is formed. Another significant contributor to ΔΔE is Val-67β, which makes contact with the bound CO (interestingly, the corresponding α-chain residue, Val-62α, is ≈1 Å further away), and is displaced on ligation. Likewise contributing to ΔΔE is Gly-64β, which is in direct contact with Val-67β. Finally, Phe-103β is an important contributor because of its contact with pyrrole 3 of the heme group [supporting information (SI) Fig. 5], which is stressed in the T structure (see below).
The more distributed character of the β-chain than the of α-chain ΔΔE values may reflect greater flexibility as a result of residue deletions (46–51 in the β-chain) in the CD loop of the α-chain. Normal mode analysis (11, 34) indicates large-amplitude motions in the CD loop. Analysis of the crystal structures has revealed greater R/T displacement of the E-helix in the β- than in the α-chains, like those exhibited by Mb (20), whose CD loop is similar to that of the Hb β-chain.
Fig. 4 shows the three regions where the main contributors to the CO affinity differential reside. When considering how the ligation energy is partitioned among different residues, it is clear that a cumulative effect results from the sum of several small terms as has been suggested by NMR studies (35). The small terms are mainly associated with the stress energy of the E- and F-helix motions. To be able to draw a quantitative picture of these results requires application of free-energy methods. These methods attempt an extensive sampling of the energetics of system (36), but this type of approach is beyond the scope of the current study.
Schematic showing the heme, proximal, and distal histidines, E- and F-helices. (Insets) Magnifications of those regions that were determined to be important for affinity modulation of CO. The arrows highlight the motions of the E- and F-helices, both are pushed downward in the deoxy state and upwards when ligation takes place.
Heme Geometry Contribution to Affinity Change.
How are the T/R energetics affected by protein-imposed distortions of the heme structure? To examine this question, we extracted the heme and its proximal ligand (SI Fig. 5) in silico and computed its QM energy from the five independent QM/MM geometries that had been optimized in the protein. Then the heme was allowed to relax by reoptimizing the geometry in the absence of protein constraints, and its energy was recalculated. As expected, relaxation energies were much larger (5–6 kcal/mol) in the T than the R structure for the singlet state of the heme, reflecting greater protein-induced strain associated with ligation in the T structure (Table 2).
Relaxation energies for the model system (kcal/mol)
To further examine the sources of geometry-derived energy changes, we divided the heme complex into subsystems associated with the CO and imidazole ligands and the heme itself, and calculated their separate contributions to the singlet-quintet energy differences (SI Table 3). The results indicate that the heme itself is the major contributor, with the imidazole ligand being the next important player, especially in the T structure (37, 38). The larger imidazole contribution to the stress comes from α-chain in the T structure (≈2 kcal/mol, twice that of the β-chain), in line with resonance Raman data showing a larger R/T reduction in the Fe-imidazole stretching frequency for the α- than for the β-chains (39).
Finally, we inspected the source of protein stresses on the heme by further subdividing it into the four pyrrole rings (see SI Fig. 5 for numbering), together with their substituents (SI Table 4). The main contributors are pyrroles 3 and 4, especially in T structure, and to a greater extent in the β- than in the α-chains. The particularly large T/R difference for the pyrrole 3 (1.3 kcal/mol) is consistent with the large contribution from Phe-103β, which packs against it, as noted above. We note also that porphyrin reconstitution studies have revealed a significant effect on ligand affinity when the vinyl substituent is replaced on pyrrole 3 but not on pyrrole 2 (40, 41). Warshel and Weiss (17) had earlier identified the importance of the vinyl group at pyrrole 3 for the computed affinity reduction in binding O2.
Conclusions
QM/MM computations provide deoxy/CO-heme energy differences within Hb that are in reasonable accord with thermodynamic data. They successfully model the reduction in ligand affinity observed for Hb between the R and T quaternary structures. About 60% of this reduction is attributable to interactions in the protein, and the remainder to strains imposed on the heme and its ligands. About half of the heme strain is associated with the porphyrin ring, and especially with pyrrole 3, which experience significant steric contacts with adjacent phenylalanine residues. Significant strain is also associated with the proximal histidine.
Consistent with experiment, the R/T energy reduction is the same for α- and β-chains, but the mechanisms are different. Although the C-terminal T constraints are indicated to be significant for both chains, proximal histidine strain is greater for the α-chains, consistent with the spectroscopic data (39). In addition, strain is focused in the FG corner of the α-chains, at the site of the important hinge quaternary contact, which is known to be critical for T state stabilization and cooperative ligand binding. These results are strongly supported by the cryogenic crystallography of Shibayama and coworkers (23). They showed that photodissociation of CO in the T structure (produced in Fe,Ni hybrids) leads to significant movement of the α-chain F-helix, including the proximal histidine that it contains, which is propagated to the FG corner. The corresponding movement in the β-chain is much smaller. On the other hand, the crystallography shows a significant motion of the β-chain Val-67β, in contact with the bound CO (23), consistent with the significant R/T energy difference computed for this residue. However, the energy differences are more broadly distributed in the β-chains, suggesting greater flexibility, which may be associated with a longer and more mobile CD loop.
Long ago Perutz proposed that T state constraints are expressed differently for α- and β-chains: proximal strain in the former and distal crowding in the latter (2). Subsequent studies have strengthened and fleshed out this hypothesis, and the present computation lends it quantitative support. It seems remarkable that such different mechanisms produce quantitatively similar reductions of ligation affinity. However, as others have pointed out (26), similar affinities in both quaternary structures are a requirement for maximum cooperativity in ligand binding and efficient oxygen delivery. Consequently there would have been evolutionary pressure for both chains to produce equivalent affinities despite their differing structures.
Methods
Initial Preparation.
Protein Data Bank (PDB) structures 1A3N (deoxy, T state) and 1BBB (carbonmonoxy, R2 state) were protonated and missing residues were filled in with the PLOP code (42). Distal histidines were assumed to be protonated at the ϵ-position, in both T and R structures (28). The protonation state of the rest of the histidines was assessed by visual inspection. The whole constrained protein was solvated in a 66 × 72 × 81 Å box of water (SPC, simple point charge model) and dynamics performed for a short 5-ps trajectory to allow the water to relax.
Ligand Migration.
Both T and R structures were allowed to relax with Protein Energy Landscape Exploration (PELE). PELE combines protein structure prediction algorithms, a ligand perturbation, and a Metropolis acceptance scheme to explore the protein–ligand energy landscape. The PELE method has been tested with applications to Mb conformational changes and ligand migration (20, 43), showing excellent conformational sampling capabilities. In Hb's T structure, placement of the CO molecule into a heme pocket was carried out by initially placing the CO next to the CD loop region, a region that we had already observed to be a preferred exit path for the ligand (43, 44), and then directing it through a series of intermediate target points until it reached the distal heme pocket. From these sampling runs, five distinct snapshots were selected as initial conformations for the QM/MM simulations. Initial conformations were generated for both α- and β-chains and the T and R structures, so that a total of 20 different structures were generated in this way.
QM/MM Simulations.
For the QM/MM simulation the system was reduced by deleting all waters beyond 10 Å from the chain of interest. The system was partitioned into two regions: the quantum (QM) region that comprised CO, heme, and proximal histidine of one chain, and the molecular mechanics (MM) region that included the rest of the protein (including the other three hemes) and water. Regions where the QM and MM region are covalently attached (residues flanking the proximal histidine) use a set of frozen orbitals specifically parametrized for the protein backbone (45).
Simulations were performed by using QSite (46) from the Schrödinger suite of programs. This protocol ensures the computation of an affinity snapshot for a given quaternary state. Quantum calculations were performed with the unrestricted density functional level of theory by using the B3LY functional and LACVP*/6–31G* basis functions. The classical (MM) region used the parameters from the OPLS 2001 force field (47).
QM Simulations.
QM calculations for the isolated model system (CO, heme, and imidazole) were performed with Jaguar (48) and the same setup as the QM region of the QM/MM simulations. Geometries of the model systems were extracted from the previously relaxed QM/MM structures described earlier. In these reduced systems the propionate groups were substituted by capping hydrogens to prevent fortuitous hydrogen bond formation with the imidazole ring. The effect of the substitution was found to be small when taking energy differences. A single-point calculation was performed to evaluate the energy of the initial geometry. Afterward, the complete model system was minimized to allow calculation of a full relaxation energy. In addition to this, the initial geometry was divided into three subsystems. Each subsystem was minimized independently by employing a constraining scheme. A constraining scheme allowed only the CO, heme, or imidazole subsystem to move during the minimization, whereas the remaining subsystems maintained their initial geometries. Finally, a relaxation energy was calculated for each subsystem (CO, heme, and imidazole) by taking the energy difference between initial and resulting geometries. The sum of the individual subsystem relaxation energies compared well with the full-relaxation energy and gave support to the assumption of independent subsystems. The same strategy was used to partition the heme into four sections containing each of the pyrroles. Independence of the subsystems was also found to be valid by comparing the full-relaxation energy of the heme with that of the combined pyrroles.
Acknowledgments
We thank Richard A. Friesner, Gary Ackers, Jiali Gao, and Arieh Warshel for their helpful discussions. Computational resources were provided by the Barcelona Supercomputing Center. This work was supported by startup funds from Washington University in St. Louis (to V.G.), the Barcelona Supercomputer Center (to V.G.), National Institutes of Health Grant GM 33576 (to T.G.S.), and Spanish MEC Grant CTQ2006-10262 (to V.G.). V.G. is an ICREA Research Professor of Life Sciences at the Barcelona Supercomputing Center.
Footnotes
- ¶To whom correspondence should be addressed. E-mail: victor.guallar{at}bsc.es
-
Author contributions: R.E.A., C.X., T.G.S., and V.G. designed research; R.E.A. and C.X. performed research; R.E.A., C.X., T.G.S., and V.G. analyzed data; R.E.A., T.G.S., and V.G. wrote the paper.
-
The authors declare no conflict of interest.
-
This article is a PNAS Direct Submission.
-
This article contains supporting information online at www.pnas.org/cgi/content/full/0706026104/DC1.
- Abbreviations:
- Mb,
- myoglobin;
- QM/MM,
- quantum mechanics/molecular mechanics;
- R,
- relaxed;
- T,
- tense.
- © 2007 by The National Academy of Sciences of the USA