Search tips
Search criteria 


Logo of scirepAboutEditorial BoardFor AuthorsScientific Reports
Sci Rep. 2017; 7: 10926.
Published online 2017 September 7. doi:  10.1038/s41598-017-11259-0
PMCID: PMC5589765

Tertiary and quaternary structural basis of oxygen affinity in human hemoglobin as revealed by multiscale simulations


Human hemoglobin (Hb) is a benchmark protein of structural biology that shaped our view of allosterism over 60 years ago, with the introduction of the MWC model based on Perutz structures of the oxy(R) and deoxy(T) states and the more recent Tertiary Two-State model that proposed the existence of individual subunit states -“r” and “t”-, whose structure is yet unknown. Cooperative oxygen binding is essential for Hb function, and despite decades of research there are still open questions related to how tertiary and quaternary changes regulate oxygen affinity. In the present work, we have determined the free energy profiles of oxygen migration and for HisE7 gate opening, with QM/MM calculations of the oxygen binding energy in order to address the influence of tertiary differences in the control of oxygen affinity. Our results show that in the α subunit the low to high affinity transition is achieved by a proximal effect that mostly affects oxygen dissociation and is the driving force of the allosteric transition, while in the β subunit the affinity change results from a complex interplay of proximal and distal effects, including an increase in the HE7 gate opening, that as shown by free energy profiles promotes oxygen uptake.


Hemoglobin (Hb) and the simpler related protein, myoglobin (Mb), are among the most thoroughly studied heme proteins, being benchmarks cases for many experimental and theoretical developments in a wide range of fields covering structural biology, bioinorganic chemistry and genetics, among others. Hb is a tetrameric protein composed of four globular (Mb-like) subunits, two of each kind, named α and β. Each subunit is composed of 8 α-helices, named A to H, connected by short non helical regions and holds the prosthetic heme group (see Fig. 1) coordinated to the protein through histidine F8 (HF8), called the proximal histidine. The other key residue related to its function is the histidine present in position E7, the distal histidine (HE7)1, 2. Hb function is tightly related to heme iron reactivity towards small molecules, specially oxygen, and even though its history spans almost a century, several questions remain open3, 4. Since Hb is a tetrameric protein, the key to its physiological role concerns the cooperatively regulated differential oxygen affinity between the T(tense, low affinity) and R (relaxed, high affinity states), and how it relates to both quaternary and tertiary structural changes.

Figure 1
Hemoglobin tetrameric alpha helical structure on the left. Schematic representation of Hb heme environment, highlighting both proximal and distal histidines, on the right.

Allosteric and structural based models used to understand oxygen affinity in Hb were born with the initial Hb crystal structures obtained by Perutz for both oxy(R) and deoxy(T) Hb states5, 6 and later developed into general allosteric MWC model by Monod et al.7. This symmetry-based model relies on the existence of an equilibrium between the two quaternary R and T states, and assumes equivalence of all four subunits. The fact that Hb contains two chemically different subunits α and β and presents a single two-fold axis of symmetry, lead to extensions of the MWC model, as well as alternative formulations such as the sequential model, introduced by Pauling8 later elaborated by Nemerthy et al. and usually known as NFK model9. Perutz also spent several years looking for structural reasons underlying the different states oxygen affinity and proposing that oxygen binding to the T state results in an in-plane movement of the heme iron, with the concomitant motion of the proximal HisF8 and the associated F-helix, resulting in the breaking of several key intersubunit salt bridges that would promote the transition to the R-state, leading to the so called stereochemical model5, 6, 10, 11. Decades later, Szabo and Karplus revisited this model and included statistical thermodynamic and energetic considerations into the stereochemical model, proposing that low oxygen affinity of the subunits in the T state, arose from the strain at the ‘allosteric core’, mainly the steric repulsion between HisF8 and the porphyrin12, 13. All these models focused on quaternary transitions and paid little attention to differences between type of subunits and tertiary transitions.

The latest -more recent- model, which explicitly includes tertiary differences, is known as the Tertiary Two State (TTS) model, as originally proposed by Henry et al.14 and was motivated by observations of the oxygen affinity in quaternary trapped states15, 16 and spectroscopic evidence of pure tertiary transitions17. The TTS model postulates that high and low affinity conformations of individual subunits, which are called “r” and “t”, exist in equilibrium within each quaternary structure, and allows -in opposition to MWC- incomplete coupling between tertiary and quaternary transitions. The TTS model recently received strong support from CO rebinding kinetic experiments in gel encapsulated Hb by Viappiani et al.18.

At this point it is important to note that although the described Hb models allow for the interpretation and assignment of oxygen binding data, and relate them with some structural features, a clear microscopic picture of how conformation of each tertiary/quaternary state control oxygen differential affinity is still missing. In this context, several basic questions like “what is the structural origin of the difference in affinity for oxygen binding to t and r?” or, “what are the relative contributions of the proximal strain and intersubunit contacts?” remain unanswered.

Chemically speaking, oxygen -and any other ligand- affinity in heme proteins, can be rationalized in terms of two processes: oxygen association, characterized by the bimolecular rate constant termed “kon”, and oxygen release, which is characterized by the dissociation rate constant “koff19. Kinetic measurements using hybrid Hbs, in which alternatively α or β hemes are replaced by porphyrin containing metals that do not bind O2 or CO (CrII, MnII, NiII and MgII), allowed estimation of each subunit quaternary dependent rate constants, yielding values of kon(α) = 11μM−1s−1, kon(β) = 5μMs−1, koff(α) = 3700 s−1 and koff(β) = 1800s−1 of the T state, and kon(α) = 36μM−1s−1, kon(β) = 76μM−1s−1, koff(α) = 16 s−1 and koff(β) = 32 s−1 for the R state20. In globins, ligand association is controlled by migration of the ligand from the bulk solvent into the active site through the protein gates and/or tunnels, a process -that as shown in our previous works, can be characterized by computing the corresponding free energy profile (FEP)2123. In particular, for Mb, we showed that the depth of the free energy well in the distal pocket, which is directly related to kon, is controlled by the HE7 “gate” conformation24. Dissociation, on the other hand, is controlled by disruption of the protein bound ligand interactions, which, as studied in previous works2527, can be accurately determined using hybrid quantum mechanics/molecular mechanics (QM/MM) based methods and directly correlated with observed koff values25. More important, both calculations, allow a deep understanding of the structural reasons that control each process and thus ultimately determine globin ligand affinity28.

In the present work we have combined the above mentioned type of classical molecular dynamics (MD) and multiscale QM/MM simulations to address the influence of tertiary structures in the control oxygen affinity, and the molecular basis of the difference between each subunit type (α vs β) and quaternary state (R vs T) oxygen affinity. Our results, discussed in the context of Hb allosteric models, show how proximal effects that govern α chain differential affinity are the driving force of the allosteric transition, and how proximal and distal effects are combined to modulate β chain affinity.

Computational methods

Starting Structures

Two representative high resolution X-ray crystallographic structures were taken as starting points for all simulations, corresponding to the deoxy (T) and oxy (R) forms of Hb, (PDB codes 2dn1 and 2dn2 respectively)29. These, as well as other lower resolution crystal structures corresponding to the same quaternary states, have been successfully used in previous computational studies to represent each quaternary structure2, 3, 3032. Standard protonation states at physiological pH were assigned to all ionizable residues (Asp, Glu, Lys and Arg) paying special attention to the protonation of histidines, which were assigned on the basis of the hydrogen-bond pattern with neighboring residues. In particular, β146 histidine is known to be a key factor in the quaternary T state stabilization, was simulated as the protonated (HIP) tautomer30. For the proximal heme-bound histidine (HF8), protonation was chosen to be in the Nδ position, since this is the protonation state that allows coordination to the iron. For distal histidine, (HE7) all three protonation states (HID, HIP and HIE) were simulated and analyzed (see below).

To build the proximal/distal single site mutants, we used the same protocol as in Capece et al.26. Briefly, HE7G was built by simply removing the HE7 side chain and adding a hydrogen atom in place of the Cα-Cβ bond. HF8Gly + Im mutant was built removing the His Cβ methylene and completing the resulting Gly and Imidazol with the corresponding hydrogen atoms. Both mutants were built starting from the corresponding oxygenated classically optimized structures for each subunit/state.

Classical Molecular Dynamics Simulations

All Molecular Dynamics simulations were performed using the PMEMD module of the Amber14 package. The Amber99SB force field33 was used for all residues but the heme, whose parameters were developed and thoroughly tested by our group in previous works27, 3436 and are available in the SI. It is important to note that our developed heme parameters are coordination and oxidation state specific and are specially designed to reproduce the experimentally observed different heme conformations and interactions with the protein matrix, allowing the performance of stable MD simulations of both R and T states of Hb (see below)28, 37. Complete starting structures were finally immersed in an octahedral box of TIP3P waters38. All simulations were performed using periodic boundary conditions and Ewald sums to treat long range electrostatic interactions39, the SHAKE algorithm to keep bonds involving hydrogen atoms at their equilibrium length40, a 2 fs time step for the integration of Newton’s equations, and the Berendsen thermostat and barostat to control the system temperature and pressure respectively41. Equilibration consisted of an energy minimization of the initial structures, followed by a slow, 2 ns long heating up to 300 K (in the NVT ensemble), and a final 300 ps long density equilibration procedure (in the NPT ensemble). Production MD runs consisted of 100 ns long trajectories. Frames were collected at 20 ps intervals, which were subsequently used to analyze the production trajectories.

QM–MM calculations

The initial structures for the QM–MM calculations corresponded to classically optimized starting structures described above, thus they remain close to the crystallographic conformation. Full hybrid QM–MM geometry optimizations were performed using a conjugate gradient algorithm. Only residues located less than 10 Å apart from the heme reactive center were allowed to move freely in the QM–MM runs. All QM–MM computations were performed at the DFT level using the generalized gradient approximation to the exchange-correlation energy proposed by Perdew et al.42 with the SIESTA code43 which has been extensively used by our group to study heme proteins. Basis functions consist of localized (numerical) pseudo-atomic orbitals, projected on a real space grid to compute the Hartree potential and exchange correlation potential matrix elements. For all atoms basis sets of double zeta plus polarization quality were employed with a pseudoatomic orbital energy shift of 25 meV and a grid cutoff of 150 Ry. Further details can be found elsewhere4446. In all cases, the iron porphyrinate, the O2 ligand and the imidazole group of the proximal histidine of one Hb subunit were selected as the quantum subsystem. The rest of the protein, together with water molecules, were treated classically.

To analyze structural effects on oxygen dissociation, we computed for each case the O2 binding energy (ΔEO2), which is defined as:

ΔEO2EHemeO2 − (EO2EHeme)

where EHeme–O2 is the energy of the oxy-form of the protein, EHeme is the energy of the deoxy-form, and EO2 is the energy of the isolated oxygen molecule. This strategy was successfully used in our previous studies, particularly for determining distal and proximal effects on the oxygen affinity in globins26. Moreover, recent works from our group showed that QMMM calculated ΔEO2 can be used to estimate koff based on an statistical model, where the ΔEO2 and koff are expressed relatively to those of a free heme.47. The corresponding results show that the accuracy of the predicted dissociation rate values is around one order of magnitude and, most important, to correctly assign the order of the predictive experimental constants.25, 27 This correlation was already applied to the study of many globins oxygen affinity25. In other words, larger ΔEO2 will definitely result in smaller koff values and vice-versa.

Free Energy profiles

In order to evaluate the oxygen association process, we calculated first the free energy profiles (FEP) for O2 migration from the solvent to the distal pocket using multiple steered molecular dynamics (MSMD), as performed by Boechi et al. for Mb24. Briefly, in MSMD an external force, described as an harmonic potential, is added to the system Hamiltonian whose equilibrium value moves at an arbitrary velocity and thus steers the system towards the desired state. Integrating the external force allows determination of the irreversible work profile W(i) along the selected reaction coordinate. Finally, FEP is obtained by exponential averaging multiple work profiles, starting from different initial equilibrated microconfigurations, employing Jarzynski’s equation.48, 49


where the brackets indicate an infinite average over the collected works (W i), k B the Boltzmann constant and T the temperature, and ΔG the free energy change. Error estimation was performed according to Bustamante et al.50, considering the FEP statistical uncertainty as the square root of the Mean Square Error (MSE), computed according to MSEσ2(N) + B2(N), where σ 2 corresponds to the variance of the calculated work profiles, B 2 to the Bias of the Jarzynski’s estimator, and N the number of MSMD trajectories. (see reference50 for details). For each profile, we performed 40 trajectories using as reaction coordinate the Fe-O distance and a 150 kcalmol−1Å−2 force constant. The pulling speed was chosen to be 0.25 Åns−1 as optimized for these type of calculations in previous works23. Due to the large extension of the reaction coordinate, simulations were divided in two regions: the first region ranges from the solvent to the first minimum, which corresponds to a secondary docking site in the migration pathway. The second region covers the reaction coordinate from the docking site until the oxygen reaches the sixth coordination position. This strategy allows to keep the ligand close to the migration pathway, avoiding unproductive trajectories. In the four cases, the resulting profiles were combined in order to overlap as best as possible the free energy minimum (corresponding to the docking site) and the slope corresponding to the barrier from the docking site to the iron bound state. To qualitatively describe the Fe-ligand bond formation when the ligands gets too close to the heme we added a Morse potential for the Fe-O bond. We used 1.8 Å as the equilibrium distance, a width control parameter of 5.6 Å−1 and a well depth of 10 kcalmol−1. It is also important to note, that although charge distribution of the heme and its ligand are altered during the binding process, and ligand parameters have been proposed and implemented51, 52 to account for this fact, our purpose is to evaluate the migration free energy profile at distances which are significantly larger than Fe-O bond where these charge redistribution effects are not expected to be of relevance, as previously shown. Thus we used unbound heme and oxygen parameters during the whole SMD process. Finally, in order to avoid unproductive trajectories due to protein rotation and free movement of the ligand in water, all the SMD simulations were started with the oxygen molecule completely solvated but next to the protein surface (at the entrance of the tunnel).

Secondly, and in order to complete the free energy landscape for O2 uptake, we also computed the FEP associated to the opening of the HE7 gate using umbrella sampling simulations, following a similar protocol as in our previous work for Mb53. The chosen reaction coordinate was the dihedral angle C-Cα-Cβ-Cγ of HE7, and was sampled from 60 degrees (closed conformation) to 160 degrees approximately (open conformation). The reaction coordinate was explored in 10 consecutive windows.


The results are organized as follows: in the first section we present QM/MM analysis of oxygenated structures and their O2 binding energy for wild type as well as proximal and distal histidine mutants. Secondly, we show the results obtained in equilibrium molecular dynamics simulations on the R and T states for the tetrameric structures. Lastly, we analyze the process of O2 entrance to the distal cavity by means of free energy calculations.

QM/MM oxygen binding energy in α and β subunits of R and T hemoglobin

Wild type (wt) Crystal and solution structures

To start unraveling the tertiary structural determinants of the differential oxygen affinity in α and β subunits of R and T Hb (from now on called αR, αT βR and βT), we computed the optimized oxygenated structure, and the oxygen binding energy (ΔEO2) for each monomer in the following crystal structures 2dn1(R) and 2dn2(T), which are commonly used as R and T quaternary references structures. ΔEO2 is related with the experimental observed koff value, and is able to accurately describe and dissect the energy required to break the bound oxygen-protein (both Fe-O2 and intermolecular protein O2 interactions). Thus, proteins (or subunits) with large negative ΔEO2 display small values of koff, while those displaying smaller ΔEO2 values exhibit large koff values, as previously shown by our group for a large set of globins27. The results are presented in Table 1.

Table 1
Energetic and geometric characterization of α and β HbA subunits from both R and T crystallographic structures.

The obtained results show that, as expected, in the R state both subunits display larger ΔEO2, and thus slower ligand release rates, which are expected to result in higher affinity. Interestingly, the increase is slightly higher for the β subunits, but both subunits show similar binding energy in both tertiary states. Analysis of the structural parameters also reveal the origin, and subtle differences, in the O2 stabilization mechanisms. In the four cases, bound oxygen is forming an H-bond with the imidazole group of distal HE7, mainly as a consequence of the negative charge (qO2 < 0) acquired by the ligand upon binding due to the well known π back bonding effect, and reflected in the dHE7H-O1/O2 distances, which range between 1.8 and 2.5 Å. β subunits in both R and T states show larger charge and longer Fe-O distance, a fact that could reflect stronger distal interactions.

Concerning proximal effects, previous work from our group26 showed that the key parameter is the Fe-HisF8 distance, which corresponds to 2.12 Å for a geometry optimized isolated heme with the same method. Longer distances indicate that the protein is pushing the histidine away, reducing the O2 affinity, while the opposite effect is observed for shorter distances. Subunit αT shows the largest Fe-His distance (2.19 Å), which indicates that in this case structure is pulling HF8 away (possibly reducing oxygen affinity). In the αR state HF8 is allowed to come closer, reaching its optimal distance. In contrast, no difference is observed between the βR and βT structures, the histidine being slightly closer than in the isolated active site in both states. Finally, analysis of HF8 rotational angle and/or HF8Nδ H-bonds, two factors which have also shown to further modulate oxygen affinity, showed no significant differences between different tertiary structures.

Proximal HF8 and Distal HE7 contributions

To further evaluate how proximal or distal effects contribute to oxygen affinity in each subunit/state we computed ΔEO2 for both HE7G and HF8G + Im mutants (see methods), a strategy which has already proven successful to identify and quantify them2628. In the case of HE7, the mutation simply eliminates the distal H-bond with the bound ligand, while in the case of HF8G + Im, the added free imidazol allows better comparison with wt protein, but destroys protein capacity to control its effect. The results are presented in Table 2.

Table 2
Oxygen affinity, selected geometrical parameters and charges for HE7G and and HF8G + Im mutants in α y β subunits, for the T and R states.

Beginning the analysis with the α subunits, it is clear that distal HE7 plays almost no role in stabilizing the bound oxygen. The proximal effect, however, as revealed by the HF8G + Im mutant is striking. While in the αR state, the mutant has a 3 kcal/mol lower affinity, consistent with the close HF8-Fe distance (push effect), in αT state mutant increases the oxygen binding energy in ca. 4 kcal/mol, corroborating that in this case tertiary structure is lowering the affinity by pulling HF8 away. Thus, it is clear that the decrease in the oxygen release rate in the T-to-R transition in the α subunits is direct and exclusively related to the proximal HF8. For the β subunits the picture is more complex since both states show distal and proximal contributions to an increase in affinity, as revealed by the decrease in binding energy in all cases. The decrease is however, slightly larger for the distal effect in the R state (2 vs 1 kcal/mol) partially explaining the overall increase.

Dynamical behavior of both Hb quaternary states

To complete the above described tertiary state picture, and analyze dynamics of Hb in both quaternary states (R and T) we performed 100 ns classical MD simulations of both tetramers in the oxygenated as well as their free forms. All four systems are stable during the simulation timescale and remain close to the starting crystal structure (RSMD < 1.8 Å). Analysis of Hb mobility using RMSF plots, shows that α subunits present significantly lower mobility than β subunits, in both R and T states (see figure S1). Distal HE7 hydrogen bond analysis shows that for the T states, in the β subunit the hydrogen bond is slightly more flexible (and possibly weaker), while for the R state no significant differences are observed (see Fig. S2 and Fig. S3).

Oxygen uptake

Free energy profiles of oxygen uptake

We now turn our attention to the oxygen uptake process, experimentally characterized by the kinetic association constant (kon). As shown in our previous work for Mb, oxygen uptake is tightly connected to the distal HE7 conformation and protonation state. In Mb, HE7 has been shown to be able to adopt both a so-called closed and open conformations, which show a marked difference in the free energy associated to oxygen uptake24. To analyze this effect we computed the free energy profile (FEP) of oxygen entry in each subunit/state in both HE7 open and closed conformations.

The results presented in Fig. 2 show that, as for Mb, in all four Hb analyzed cases, there is a clear difference between the FEP in the closed and open states, the former displaying barriers in the 4-6 kcal/mol range, while the latter shows very small or no barrier at all. Moreover, also as observed in Mb, HE7 open state results in the creation of a small hydrophobic cavity that draws oxygen inside Hb active site with a a −2 kcal/mol decrease in the free energy. No significant differences are observed between the subunits/states (see S4 and S5 for work distributions and convergence analysis).

Figure 2
On top, schematic view of the open and closed states of the HE7 gate. On the bottom, free energy profiles for O2 migration from the solvent to the distal site in the closed state (left), with HE7 in the HIE protonation state, and in the open state (right), ...

Free energy profiles of HE7 opening

Given that, as shown above, oxygen uptake is dependent on HE7 opening, we decided to analyze the free energy associated to this process in each Hb subunit/state, and in all cases considering all three possible histidine tautomers, HID (protonated in δ nitrogen), HIE (protonated in ε nitrogen), and HIP (doubly protonated, with a net positive charge). The resulting FEPs, presented in Fig. 3, show very interesting trends. It is important to make explicit that each free energy profile is independent from the others, and we have translated all closed conformations to the same zero. Those FEP allow us to compare the open and closed conformations, and the barrier between them, and not the relative free energies among all subunits and distal histidine protonation states.

Figure 3
Free energy profiles associated to the HE7 gate opening process in the R (left) and T (right) allosteric states.

In the T state, both subunits show a similar behavior, with HIE being preferentially closed and with high opening barrier, HID displaying a moderate barrier, and finally the protonated HIP state showing a preference for the open state and smaller barrier. Also important, is that in the T state for all three tautomers β subunits show a favored open state compared to α subunits. In the R state the picture is slightly more complex, α subunits are again predominantly closed in HIE/HID tautomers, and even for the protonated HE7, with a slight preference for the open state, the barrier is quite high. In contrast the β subunits show in all cases small (or completely lack of) opening barrier and the open state is the preferred one in HID and HIP, displaying same free energy in HID.

In summary, while the T-to-R transition does not significantly affect open-close equilibrium in the α subunits, the open fraction -and therefore oxygen uptake- are significantly increased in the β subunits in all considered cases.


The relevance of human Hb as the paradigmatic, allosteric structurally characterized protein, lead to the performance of the first computational studies of the subject almost 40 years ago by recent Nobel laureates Karplus54 and Warshel55. Since those times, many Hb simulations were performed, mostly focusing on the allosteric process30, 56, paying less attention to how oxygen affinity is differentially regulated. Our results, together with previous works from the group on different members of the globin family, now allow a deeper understanding of the tertiary factors underlying each subunit type and state control of oxygen affinity. A qualitative summary of our results is presented in Fig. 4. Starting from T state, it is evident that in the α subunit low affinity is achieved by a proximal effect that significantly promotes a destabilization of bound oxygen leading to a high koff as experimentally observed. For the β subunit the low affinity results from a complex interplay of proximal, distal and intrinsic effects, which also leads to a high dissociation rate. Also, in the T state HisE7 open conformation is not favored, and the high barrier for oxygen entry in the closed conformation results in slow oxygen entry. The quaternary T-to-R transition, is translated into several effects that altogether increase oxygen affinity. In the α subunit, as expected, the largest and almost exclusive change is related to the proximal HF8 position which now is properly positioned to significantly increase oxygen binding energy. In the β subunit, both proximal and distal effects contribute to the increase in oxygen binding energy consistent with the observed lower koff, but also HE7 significantly increases the population of the open state as evidenced by the change in the close-to-open free energy profiles, which results in a net increase in oxygen uptake, again in agreement with the observed increase in kon.

Figure 4
Tertiary mechanisms controlling oxygen affinity in human Hb monomers. Red arrows indicate a negative effect and blue arrows, positive effect. Proximal effects are indicated in the HF8 and distal effects, in the HE7. Curved green arrows in the HE7 show ...

Our results can also be analyzed in the context of the early proposal from Perutz57 and more recent Resonance Raman studies of hybrid Hbs58, computational studies of the allosteric transition32, QM/MM studies of CO binding31 and the already mentioned TTS model14 which together, allowed to establish a connection between Hb reactivity in the context of tertiary and quaternary structural changes. Taken together, those results showed that heme reactivity is controlled exclusively by tertiary forces defining the r and t states, and that quaternary T and R conformations must work by biasing tertiary populations. They also showed that low- and high- affinity tertiary structures, presumably t and r, can coexist within each quaternary state. Our results, show that from an oxygen affinity perspective, crystal structures of each subunit in each quaternary state (αR, αT, βR and βT) are good representatives of the required tertiary low (t) and high (r) affinity structures. Moreover, the changes in affinity are predominantly due to changes in the heme and proximal HF8 -specially in the α subunits- changes that are expected to occur even if quaternary structure is frozen, thus allowing high affinity r-like state to be able to exist in a quaternary T state framework, and vice-versa. Our results, show that from an oxygen affinity perspective, the presently studied crystal structures of each subunit in each quaternary state (αR, αT, βR and βT) are good representatives of the required tertiary low (t) and high (r) affinity structures. In other words, since we observed (and rationalized) known affinity differences between α/βR and α/βT we can reasonably assign them to the “t” and “r” states of the TTS model. And since the changes in affinity are predominantly due to relatively small structural rearrangements in the heme and proximal HF8 (specially in the α subunits). These changes are expected to occur even if quaternary structure is frozen, thus allowing high affinity r-like state to be able to exist in a quaternary T state framework, and vice-versa.

Concerning the observed differences between α and β subunits, early ideas from Perutz, already pointed to proximal effects in the α chains, and distal effects in the β chains, as main contributors to the T-to-R change in affinity. These ideas, received further support from mutagenesis studies of HE7 and VE1159, 60, crystal structures of NO bound Hb which highlighted the proximal strain in αT61, 62 and above mentioned Raman Studies and QM/MM studies58. It is also interesting to note that, in spite of their underlying different mechanisms, both chains show a very similar change in oxygen affinity along the T-to-R transition. This chain equivalence, which was already observed in QM/MM calculations of CO binding63 and kinetic measurements20, has been suggested as a fundamental requirement for maximum cooperativity. In this context, our data suggest that the smaller proximal effect of the β chains, which allows them to relax -and switch reactivity- faster58 is compensated by the distal effect in both koff and kon, possibly once full βT to βR transition is achieved.

Also, in relation to the affinity and ligand binding kinetic experiments performed in quaternary trapped structures which supported the TTS mode18, 64, our data suggest that most of the change in each subunit oxygen affinity between r and t states, is related to proximal and heme conformation effects, which can be considered small (in terms of atomic displacements) compared to the overall tertiary and quaternary transitions. Therefore, its is not difficult to envisage that a quaternary fixed R state Hb, can allow individual subunits to relax the heme and its proximal His to a t-like state, thus displaying low affinity, as experimentally observed.

Finally, we can also analyze how differential oxygen affinity controlling factors work to achieve positive cooperativity and trigger the R-to-T conformational transition, the hallmarks of Hb function. Starting from deoxy Hb in the T state, first molecular event that takes place when Hb is transferred to a high PO2 environment is partial oxygenation. In particular, oxygenation of the α chains most likely results in αT  αR transition due to pull effect of the bound oxygen to proximal HF8 (dFe-NHis distance diminishes from 2.20 to 2.12 Å, and concomitant 0.15 Å Fe in plane displacement). The transition increases oxygen binding energy that results in significant slower release rate. The change now possibly propagates to the β chains, which through the transition, increase their HE7 -gate- opening as revealed by the corresponding free energy profiles, thus promoting oxygen uptake. Finally, oxygen binding to the β chains, relaxes its heme and proximal strains, lowering ligand release and further increasing the affinity. In the reverse scenario, when oxyHb is found in low PO2 environment, partially deoxygenation -specially in the α chains-, promotes remaining oxygenated β subunits to change to the t state, decreasing oxygen binding energy and thus promoting its release. In this context, our results not only provide an explanation of how cooperativity is achieved but also point to possible key events that should drive the T-to-R transition.


Using state of the art MD and QM/MM simulations we have determined the tertiary factors that differentially control subunit (α vs β) and allosteric state (R vs T) dependent oxygen affinity in human hemoglobin, highlighting the relevance of proximal effects in the α chains to drive the transition. Our results further support and strengthen the Tertiary Two State model, by showing the key tertiary factors defining the tertiary “r” and “t” states, and their coexistence in both quaternary states.

Electronic supplementary material


This work was supported in part by grants from the Universidad de Buenos Aires (UBACYT project 20020150200091BA), Agencia Nacional de Promoción Cientıfica y Tecnológica (project PICT 2016-0568), CONICET (PIP 2015-2017 112201 501003 03CO). M.B. holds a CONICET Ph.D. Fellowship. M.A.M., D.A.E., and L.C. are members of CONICET. The authors thank the Centro de Computación de Alto Rendimiento (CeCAR, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires) for providing computational resources.

Author Contributions

Author Contributions

M.B. and A.A.P. performed the simulations. M.B. and L.C. prepared all figures. M.M. and L.C. wrote the main manuscript. All authors reviewed the manuscript.


Competing Interests

The authors declare that they have no competing interests.


Electronic supplementary material

Supplementary information accompanies this paper at doi:10.1038/s41598-017-11259-0

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Contributor Information

Luciana Capece,

Marcelo A. Martí, moc.liamg@olecram.itram.


1. Lukin JA, Ho C. The Structure−Function Relationship of Hemoglobin in Solution at Atomic Resolution. Chem. Rev. 2004;104:1219–1230. doi: 10.1021/cr940325w. [PubMed] [Cross Ref]
2. Yuan Y, Tam MF, Simplaceanu V, Ho C. A New Look at Hemoglobin Allostery. Chem. Rev. 2015;115:1702–1724. doi: 10.1021/cr500495x. [PMC free article] [PubMed] [Cross Ref]
3. Eaton WA, Henry ER, Hofrichter J, Mozzarelli A. Is cooperative oxygen binding by hemoglobin really understood? Nat. Struct. Mol. Biol. 1999;6:351–358. doi: 10.1038/7586. [PubMed] [Cross Ref]
4. Eaton WA, et al. Evolution of allosteric models for hemoglobin. IUBMB Life. 2007;59:586–599. doi: 10.1080/15216540701272380. [PubMed] [Cross Ref]
5. Muirhead H, Perutz MF. Structure Of Hæemoglobin: A Three-Dimensional Fourier Synthesis of Reduced Human Haemoglobin at 5.5 |[angst]| Resolution. Nature. 1963;199:633–638. doi: 10.1038/199633a0. [PubMed] [Cross Ref]
6. Perutz MF, Bolton W, Diamond R, Muirhead H, Watson HC. Structure of Hæmoglobin: An X-ray Examination of Reduced Horse Hæmoglobin. Nature. 1964;203:687–690. doi: 10.1038/203687a0. [PubMed] [Cross Ref]
7. Wyman, J. Linked Functions and Reciprocal Effects in Hemoglobin: A Second Look. in Advances in Protein Chemistry19, 223–286 (Elsevier, 1964). [PubMed]
8. Pauling L. The Oxygen Equilibrium of Hemoglobin and Its Structural Interpretation. Proc. Natl. Acad. Sci. 1935;21:186–191. doi: 10.1073/pnas.21.4.186. [PubMed] [Cross Ref]
9. Koshland DE, Némethy G, Filmer D. Comparison of Experimental Binding Data and Theoretical Models in Proteins Containing Subunits*. Biochemistry (Mosc.) 1966;5:365–385. doi: 10.1021/bi00865a047. [PubMed] [Cross Ref]
10. Perutz MF, Fermi G, Luisi B, Shaanan B, Liddington RC. Stereochemistry of Cooperative Mechanisms in Hemoglobin. Cold Spring Harb. Symp. Quant. Biol. 1987;52:555–565. doi: 10.1101/SQB.1987.052.01.063. [PubMed] [Cross Ref]
11. Perutz MF, Wilkinson AJ, Paoli M, Dodson GG. The stereochemical mechanism of the cooperative effects in hemoglobin revisited. Annu. Rev. Biophys. Biomol. Struct. 1998;27:1–34. doi: 10.1146/annurev.biophys.27.1.1. [PubMed] [Cross Ref]
12. Lee AW, Karplus M. Structure-specific model of hemoglobin cooperativity. Proc. Natl. Acad. Sci. 1983;80:7055–7059. doi: 10.1073/pnas.80.23.7055. [PubMed] [Cross Ref]
13. Lee AW, Karplus M, Poyart C, Bursaux E. Analysis of proton release in oxygen binding by hemoglobin: implications for the cooperative mechanism. Biochemistry (Mosc.) 1988;27:1285–1301. doi: 10.1021/bi00404a031. [PubMed] [Cross Ref]
14. Henry ER, Bettati S, Hofrichter J, Eaton WA. A tertiary two-state allosteric model for hemoglobin. Biophys. Chem. 2002;98:149–164. doi: 10.1016/S0301-4622(02)00091-1. [PubMed] [Cross Ref]
15. Mozzarelli A, Rivetti C, Rossi GL, Henry ER, Eaton WA. Crystals of haemoglobin with the T quaternary structure bind oxygen noncooperatively with no Bohr effect. Nature. 1991;351:416–419. doi: 10.1038/351416a0. [PubMed] [Cross Ref]
16. Rivetti C, Mozzarelli A, Rossi GL, Henry ER, Eaton WA. Oxygen binding by single crystals of hemoglobin. Biochemistry (Mosc.) 1993;32:2888–2906. doi: 10.1021/bi00062a021. [PubMed] [Cross Ref]
17. Lyons KB, Friedman JM, Fleury PA. Nanosecond transient Raman spectra of photolysed carboxyhaemoglobin. Nature. 1978;275:565–566. doi: 10.1038/275565a0. [PubMed] [Cross Ref]
18. Viappiani C, et al. New insights into allosteric mechanisms from trapping unstable protein conformations in silica gels. Proc. Natl. Acad. Sci. 2004;101:14414–14419. doi: 10.1073/pnas.0405987101. [PubMed] [Cross Ref]
19. Tsai A, Martin E, Berka V, Olson JS. How Do Heme-Protein Sensors Exclude Oxygen? Lessons Learned from Cytochrome c′, Nostoc puntiforme Heme Nitric Oxide/Oxygen-Binding Domain, and Soluble Guanylyl Cyclase. Antioxid. Redox Signal. 2012;17:1246–1263. doi: 10.1089/ars.2012.4564. [PMC free article] [PubMed] [Cross Ref]
20. Unzai S, Eich R, Shibayama N, Olson JS, Morimoto H. Rate Constants for O2 and CO Binding to the α and β Subunits within the R and T States of Human Hemoglobin. J. Biol. Chem. 1998;273:23150–23159. doi: 10.1074/jbc.273.36.23150. [PubMed] [Cross Ref]
21. Bidon-Chanal A, et al. Ligand-induced dynamical regulation of NO conversion in Mycobacterium tuberculosis truncated hemoglobin-N. Proteins Struct. Funct. Bioinforma. 2006;64:457–464. doi: 10.1002/prot.21004. [PubMed] [Cross Ref]
22. Boechi L, et al. Structural determinants of ligand migration in Mycobacterium tuberculosis truncated hemoglobin O. Proteins Struct. Funct. Bioinforma. 2008;73:372–379. doi: 10.1002/prot.22072. [PubMed] [Cross Ref]
23. Forti F, Boechi L, Estrin DA, Marti MA. Comparing and combining implicit ligand sampling with multiple steered molecular dynamics to study ligand migration processes in heme proteins. J. Comput. Chem. 2011;32:2219–2231. doi: 10.1002/jcc.21805. [PubMed] [Cross Ref]
24. Boechi L, et al. Hydrophobic Effect Drives Oxygen Uptake in Myoglobin via Histidine E7. J. Biol. Chem. 2013;288:6754–6762. doi: 10.1074/jbc.M112.426056. [PMC free article] [PubMed] [Cross Ref]
25. Bustamante JP, et al. Evolutionary and Functional Relationships in the Truncated Hemoglobin Family. PLOS Comput. Biol. 2016;12:e1004701. doi: 10.1371/journal.pcbi.1004701. [PMC free article] [PubMed] [Cross Ref]
26. Capece L, Marti MA, Crespo A, Doctorovich F, Estrin DA. Heme Protein Oxygen Affinity Regulation Exerted by Proximal Effects. J. Am. Chem. Soc. 2006;128:12455–12461. doi: 10.1021/ja0620033. [PubMed] [Cross Ref]
27. Marti MA, et al. Dioxygen affinity in heme proteins investigated by computer simulation. J. Inorg. Biochem. 2006;100:761–770. doi: 10.1016/j.jinorgbio.2005.12.009. [PubMed] [Cross Ref]
28. Capece L, et al. Small ligand–globin interactions: Reviewing lessons derived from computer simulation. Biochim. Biophys. Acta BBA - Proteins Proteomics. 2013;1834:1722–1738. doi: 10.1016/j.bbapap.2013.02.038. [PubMed] [Cross Ref]
29. Park S-Y, Yokoyama T, Shibayama N, Shiro Y, Tame JRH. 1.25 Å Resolution Crystal Structures of Human Haemoglobin in the Oxy, Deoxy and Carbonmonoxy Forms. J. Mol. Biol. 2006;360:690–701. doi: 10.1016/j.jmb.2006.05.036. [PubMed] [Cross Ref]
30. Hub JS, Kubitzki MB, de Groot BL. Spontaneous Quaternary and Tertiary T-R Transitions of Human Hemoglobin in Molecular Dynamics Simulation. PLoS Comput. Biol. 2010;6:e1000774. doi: 10.1371/journal.pcbi.1000774. [PMC free article] [PubMed] [Cross Ref]
31. Alcantara RE, Xu C, Spiro TG, Guallar V. A quantum-chemical picture of hemoglobin affinity. Proc. Natl. Acad. Sci. 2007;104:18451–18455. doi: 10.1073/pnas.0706026104. [PubMed] [Cross Ref]
32. Fischer S, Olsen KW, Nam K, Karplus M. Unsuspected pathway of the allosteric transition in hemoglobin. Proc. Natl. Acad. Sci. 2011;108:5608–5613. doi: 10.1073/pnas.1011995108. [PubMed] [Cross Ref]
33. Lindorff-Larsen, K. et al. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins Struct. Funct. Bioinforma. NA-NA doi:10.1002/prot.22711 (2010). [PMC free article] [PubMed]
34. Capece L, Estrin DA, Marti MA. Dynamical Characterization of the Heme NO Oxygen Binding (HNOX) Domain. Insight into Soluble Guanylate Cyclase Allosteric Transition Biochemistry (Mosc.) 2008;47:9416–9427. doi: 10.1021/bi800682k. [PubMed] [Cross Ref]
35. Bidon-Chanal, A., Martí, M. A., Estrín, D. A. & Luque, F. J. Exploring the Nitric Oxide Detoxification Mechanism of Mycobacterium tuberculosis Truncated Haemoglobin N. in SelfOrganization of Molecular Systems (eds. Russo, N., Antonchenko, V. Y. & Kryachko, E. S.) 33–47 (Springer Netherlands, 2009).
36. Nadra AD, Martí MA, Pesce A, Bolognesi M, Estrin DA. Exploring the molecular basis of heme coordination in human neuroglobin. Proteins Struct. Funct. Bioinforma. 2008;71:695–705. doi: 10.1002/prot.21814. [PubMed] [Cross Ref]
37. Ramírez CL, et al. Coarse-Grained Simulations of Heme Proteins: Validation and Study of Large Conformational Transitions. J. Chem. Theory Comput. 2016;12:3390–3397. doi: 10.1021/acs.jctc.6b00278. [PubMed] [Cross Ref]
38. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983;79:926–935. doi: 10.1063/1.445869. [Cross Ref]
39. Particle mesh Ewald: An N[center dot]log(N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092 (1993).
40. Ryckaert J-P, Ciccotti G, Berendsen HJ. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977;23:327–341. doi: 10.1016/0021-9991(77)90098-5. [Cross Ref]
41. Berendsen HJC, Postma JPM, van Gunsteren WF, DiNola A, Haak JR. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984;81:3684–3690. doi: 10.1063/1.448118. [Cross Ref]
42. Perdew JP, Burke K, Ernzerhof M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996;77:3865–3868. doi: 10.1103/PhysRevLett.77.3865. [PubMed] [Cross Ref]
43. Soler JM, et al. The SIESTA method for ab initio order- N materials simulation. J. Phys. Condens. Matter. 2002;14:2745–2779. doi: 10.1088/0953-8984/14/11/302. [Cross Ref]
44. Crespo A, et al. Theoretical Study of the Truncated Hemoglobin HbN: Exploring the Molecular Basis of the NO Detoxification Mechanism. J. Am. Chem. Soc. 2005;127:4433–4444. doi: 10.1021/ja0450004. [PubMed] [Cross Ref]
45. Martí MA, Capece L, Crespo A, Doctorovich F, Estrin DA. Nitric Oxide Interaction with Cytochrome c′ and Its Relevance to Guanylate Cyclase. Why Does the Iron Histidine Bond Break? J. Am. Chem. Soc. 2005;127:7721–7728. doi: 10.1021/ja042870c. [PubMed] [Cross Ref]
46. Martí MA, Crespo A, Bari SE, Doctorovich FA, Estrin DA. QM−MM Study of Nitrite Reduction by Nitrite Reductase of Pseudomonas aeruginosa. J. Phys. Chem. B. 2004;108:18073–18080. doi: 10.1021/jp048807r. [Cross Ref]
47. Bustamante JP, et al. A quantitative model for oxygen uptake and release in a family of hemeproteins. Bioinformatics. 2016;32:1805–1813. doi: 10.1093/bioinformatics/btw083. [PubMed] [Cross Ref]
48. Jarzynski C. Nonequilibrium Equality for Free Energy Differences. Phys. Rev. Lett. 1997;78:2690–2693. doi: 10.1103/PhysRevLett.78.2690. [Cross Ref]
49. Ramírez, C. L., Martí, M. A. & Roitberg, A. E. Steered Molecular Dynamics Methods Applied to Enzyme Mechanism and Energetics. in Methods in Enzymology578, 123–143 (Elsevier, 2016). [PubMed]
50. Gore J, Ritort F, Bustamante C. Bias and error in estimates of equilibrium free-energy differences from nonequilibrium measurements. Proc. Natl. Acad. Sci. 2003;100:12564–12569. doi: 10.1073/pnas.1635159100. [PubMed] [Cross Ref]
51. Nutt DR, Meuwly M. Theoretical Investigation of Infrared Spectra and Pocket Dynamics of Photodissociated Carbonmonoxy Myoglobin. Biophys. J. 2003;85:3612–3623. doi: 10.1016/S0006-3495(03)74779-1. [PubMed] [Cross Ref]
52. Nutt DR, Meuwly M. CO migration in native and mutant myoglobin: Atomistic simulations for the understanding of protein function. Proc. Natl. Acad. Sci. USA. 2004;101:5998–6002. doi: 10.1073/pnas.0306712101. [PubMed] [Cross Ref]
53. Arcon JP, Rosi P, Petruk AA, Marti MA, Estrin DA. Molecular Mechanism of Myoglobin Autoxidation: Insights from Computer Simulations. J. Phys. Chem. B. 2015;119:1802–1813. doi: 10.1021/jp5093948. [PubMed] [Cross Ref]
54. Gelin BR, Karplus M. Mechanism of tertiary structural change in hemoglobin. Proc. Natl. Acad. Sci. 1977;74:801–805. doi: 10.1073/pnas.74.3.801. [PubMed] [Cross Ref]
55. Warshel A. Energy-structure correlation in metalloporphyrins and the control of oxygen binding by hemoglobin. Proc. Natl. Acad. Sci. 1977;74:1789–1793. doi: 10.1073/pnas.74.5.1789. [PubMed] [Cross Ref]
56. Mouawad L, Perahia D, Robert CH, Guilbert C. New Insights into the Allosteric Mechanism of Human Hemoglobin from Molecular Dynamics Simulations. Biophys. J. 2002;82:3224–3245. doi: 10.1016/S0006-3495(02)75665-8. [PubMed] [Cross Ref]
57. Perutz MF. Stereochemistry of Cooperative Effects in Haemoglobin: Haem–Haem Interaction and the Problem of Allostery. Nature. 1970;228:726–734. doi: 10.1038/228726a0. [PubMed] [Cross Ref]
58. Jones EM, et al. Differential Control of Heme Reactivity in Alpha and Beta Subunits of Hemoglobin: A Combined Raman Spectroscopic and Computational Study. J. Am. Chem. Soc. 2014;136:10325–10339. doi: 10.1021/ja503328a. [PMC free article] [PubMed] [Cross Ref]
59. Nagai K, et al. Distal residues in the oxygen binding site of haemoglobin studied by protein engineering. Nature. 1987;329:858–860. doi: 10.1038/329858a0. [PubMed] [Cross Ref]
60. Aj M, Js O, Jp R. J, T. & K, N. The assignment of carbon monoxide association rate constants to the alpha and beta subunits in native and mutant human deoxyhemoglobin tetramers. J. Biol. Chem. 1991;266:21631–21639. [PubMed]
61. Martí MA, Scherlis DA, Doctorovich FA, Ordejón P, Estrin DA. Modulation of the NO trans effect in heme proteins: implications for the activation of soluble guanylate cyclase. JBIC J. Biol. Inorg. Chem. 2003;8:595–600. doi: 10.1007/s00775-003-0452-9. [PubMed] [Cross Ref]
62. Chan N-L, Kavanaugh JS, Rogers PH, Arnone A. Crystallographic Analysis of the Interaction of Nitric Oxide with Quaternary-T Human Hemoglobin †, ‡ Biochemistry (Mosc.) 2004;43:118–132. doi: 10.1021/bi030172j. [PubMed] [Cross Ref]
63. Lucas MF, Guallar V. An Atomistic View on Human Hemoglobin Carbon Monoxide Migration Processes. Biophys. J. 2012;102:887–896. doi: 10.1016/j.bpj.2012.01.011. [PubMed] [Cross Ref]
64. Viappiani C, et al. Experimental basis for a new allosteric model for multisubunit proteins. Proc. Natl. Acad. Sci. 2014;111:12758–12763. doi: 10.1073/pnas.1413566111. [PubMed] [Cross Ref]

Articles from Scientific Reports are provided here courtesy of Nature Publishing Group