1. Structural Biology and Molecular Biophysics
Download icon

Valid molecular dynamics simulations of human hemoglobin require a surprisingly large box size

  1. Krystel El Hage
  2. Florent Hédin
  3. Prashant K Gupta
  4. Markus Meuwly  Is a corresponding author
  5. Martin Karplus  Is a corresponding author
  1. University of Basel, Switzerland
  2. Harvard University, Cambridge, United States
  3. ISIS, Université de Strasbourg, France
Research Article
  • Cited 1
  • Views 1,137
  • Annotations
Cite this article as: eLife 2018;7:e35560 doi: 10.7554/eLife.35560

Abstract

Recent molecular dynamics (MD) simulations of human hemoglobin (Hb) give results in disagreement with experiment. Although it is known that the unliganded (T0) and liganded (R4) tetramers are stable in solution, the published MD simulations of T0 undergo a rapid quaternary transition to an R-like structure. We show that T0 is stable only when the periodic solvent box contains ten times more water molecules than the standard size for such simulations. The results suggest that such a large box is required for the hydrophobic effect, which stabilizes the T0 tetramer, to be manifested. Even in the largest box, T0 is not stable unless His146 is protonated, providing an atomistic validation of the Perutz model. The possibility that extra large boxes are required to obtain meaningful results will have to be considered in evaluating existing and future simulations of a wide range of systems.

https://doi.org/10.7554/eLife.35560.001

Introduction

Human hemoglobin is the paradigmatic model system for cooperativity in proteins. It transports oxygen from the lungs to the tissues and it is composed of two identical α-chains (α1 and α2) and two identical βchains (β1 and β2). They form two identical dimers (α1β1 and α2β2), whose relative orientation differs significantly in the unliganded (T0) and liganded (R4) tetramer. (Baldwin and Chothia, 1979) Although there is a vast literature on Hb and its cooperative mechanism (Schechter, 2008), how it functions at the atomistic level is still not fully understood (Cui and Karplus, 2008). The first insight into the mechanism was obtained from the low-resolution structures (5.5 Å) of the hemoglobin tetramer (Muirhead and Perutz, 1963), which showed that the heme groups were too distant to be able to interact directly. Monod, Wyman and Changeux formulated the allosteric (MWC) model (Monod et al., 1965; Cui and Karplus, 2008) based on the structural transition between two quaternary structures (T and R) to explain the indirect interaction between the heme groups required for cooperative oxygen binding.

Higher resolution (2.8 Å) X-ray structures of unliganded and liganded hemoglobin (Perutz, 1970) confirmed that there exist two quaternary structures (deoxy (T0) and oxy (R4)) for the tetramer and two tertiary structures for each individual chain (liganded and unliganded). Based on the structural results, as well as mutant data, Perutz, 1970) proposed a stereochemical mechanism for cooperativity, in which salt bridges (some with ionizable protons in the neutral pH range) provide the link between ligand-induced tertiary changes and the relative stability of the two quaternary structures. Shortly afterwards, the elements of the Perutz mechanism were incorporated into a structure-based statistical mechanical model (Szabo and Karplus, 1972). The model provides a quantitative framework for the effects of specific tertiary structural changes induced by ligand binding on the relative stability of the T and R structures (Szabo and Karplus, 1972; Gelin and Karplus, 1977). The shift of the equilibrium from T to R as a function of ligand concentration results in the sigmoidal (cooperative) ligand binding curve.

Several papers using different force fields and simulation conditions have been published recently (Hub et al., 2010; Yusuff et al., 2012) describing molecular dynamics (MD) simulations of the T and R states, including T0 and R4 for which 1.25 Å resolution X-ray structures are available (T0, 2DN2; R4, 2DN3) (Park et al., 2006). Although R4 was found to be stable for several hundred nanoseconds, the T0 state was not. It was found to make a transition to an R-like state in the same time period (see also Figure 1—figure supplement 2). This occurs in spite of the fact that, experimentally, the T0 state is about seven kcal/mol more stable than the R0 state (Edelstein, 1971), which is derived from the ratio of the dissociation constants of liganded and unliganded Hb of 6.7×105 (Thomas and Edelstein, 1972).

Although the rate of the T0 to R0 transition has not been measured directly, it can be estimated from the experimentally determined R0 to T0 transition rates and the R0/T0 equilibrium constant. For the unliganded hemoglobin tetramer, the R0 to T0 rate is about 20 μs (Sawicki and Gibson, 1976). With the equilibrium constant of [T0/R0] equal to 6.7×105 (Edelstein, 1971; Thomas and Edelstein, 1972), the T0 to R0 rate is estimated to be on the order of seconds with the major contribution to the activation barrier arising from the equilibrium free energy difference (7 kcal/mol) between T0 and R0. There is, thus, a striking disagreement between the transition time observed in the simulations and the estimate from experiment.

Results and discussion

The instability of T0 in the published simulations raises a fundamental question: What is wrong with them? In search for an answer, we focused on the hydrophobic effect (Chothia et al., 1976; Lesk et al., 1985), which arises from the disruption of the bulk water hydrogen bond network around nonpolar groups (Rossky et al., 1979; Cheng and Rossky, 1998). The theoretical analysis of Chandler and coworkers (Chandler, 2005; Chandler and Varilly, 2012) indicated that for large molecules, there was a ‘dewetting’ phenomenon that stabilizes a more compact structure. Chothia et al. (1976) noted that "A larger surface area is buried in deoxy- than in methemoglobin as a result of tertiary and quaternary structure changes. [..] This implies that hydrophobicity stabilizes the deoxy structure, the free energy spent in keeping the subunits in a low-affinity conformation being compensated by hydrophobic free energy due to the smaller surface area accessible to solvent.’ Such a stabilizing effect of T0 should appear naturally in a MD simulation, but evidently did not do so in the published simulations (Hub et al., 2010; Yusuff et al., 2012).

Figure 1 with 6 supplements see all
Global structural changes depending on box size.

Panel A: Overall structure of the α1β1α2β2 hemoglobin tetramer. The His146 side chains (green spheres) are specifically indicated. Panel B: Temporal change of the Cα–Cα distance between His146β1 and His146β2. The dashed lines (black, orange, red) indicate transition points for the 75, 90, and 120 box, respectively. Cyan and blue arrows indicate the values of the corresponding observable found for the deoxy T0 (2DN2) and oxy R4 (2DN3) states, respectively.

https://doi.org/10.7554/eLife.35560.002

Having exhausted other possibilities (see Appendix for more details), we wondered whether the box size used for the MD simulations might be too small for the hydrophobic stabilization to be manifested. Because most of the simulation time is consumed by the waters, rather than the protein itself, which is of primary interest, a ‘lore’ has grown up about the minimal box size that can be used with periodic boundary conditions to carry out meaningful simulations. The standard requirement is that there must be at least five water molecules between any protein atom and the box boundary. (Brooks et al., 1983, 2009) The box we first used was 75 Å; the T0 tetramer dimensions are approximately 54×49×50 Å, and there were 10,763 water molecules, as well as enough Na+ and Cl ions to yield a 0.15 m/L molar concentration. All MD simulations were done in the NPT ensemble (see SI for details).

To investigate the possibility that larger boxes were required for stabilizing T0, we carried out 1μs simulations with four cubic boxes, 75 Å (10,543 water molecules) 90 Å (20,756 water molecules), 120 Å (53,287) and 150 Å (105,073), see Figure 1—figure supplement 1. In all these simulations, His146β1 and His146β2, which play an essential role in the Perutz model (Perutz, 1970), were protonated. A comprehensive overview of the structural changes observed in the simulations of the four box sizes is provided in Figure 1. The essential result is that T0 is stable for the entire simulation in the 150 Å box, while it is not in the smaller boxes (for details, see SI). Figure 2 shows the structures obtained at the end of the 1 μs simulations, superposed on the X-ray structure that is more similar; that is, the oxy R4 structure (2DN3) for the 90 and 120 Å simulations and the deoxy T0 structure (2DN2) for the 150 Å box simulation (see also Table 1).

Global conformational rearrangement of tetrameric hemoglobin as a function of the box size.

(A) Superposition of the 2DN3 structure and the Hb structure in the 90 Å box at t=1000 ns. (B) Superposition of the 2DN3 structure and the Hb structure in the 120 Å box at t=1000 ns. (C) Superposition of the 2DN2 structure and the structure in the 150 Å box at t=1000 ns. In each case the RMSDs (in Å) to the 2DN2 and 2DN3 are also given.

https://doi.org/10.7554/eLife.35560.009
Table 1
Cα RMSD (in Å) relative to the 2DN2 (T0) and 2DN3 (R4) X-ray structure (Park et al., 2006) of the end point structures at 1 μs.

In bold structure to which a computed Hb structure is closest.

https://doi.org/10.7554/eLife.35560.010
StructureCα-Cα RMSD to 2DN2Cα-Cα RMSD to 2DN3
2DN2

2.43
2DN32.43

75 Å box2.372.59
90 Å box3.350.64
120 Å box2.440.30
150 Å box0.393.09
Figure 3 with 1 supplement see all
The average number of H-bonds per water molecule <NHbonds>/molecule, from analyzing water-water hydrogen bonds, during 1 μs MD simulations for all four box sizes and for three different donor-acceptor distance cutoffs: 2.8 (strong, top panels), 3.0 (medium, middle panels) and 3.3 Å (weak, bottom panels).

The solid black lines are averages for time intervals corresponding to the lifetime of each state in each of the simulation boxes (see Figure 1): For the 75 Å box: 0 to 140 ns (first transition), 140 to 530 ns (second transition) and 530 to 1000 ns. For the 90 Å box: 0 to 470 ns, 470 to 780 ns and 780 to 1000 ns. For the 120 Å box: 0 to 620 ns, 620 to 920 ns, and 920 to 1000 ns. For the 150 Å box, the average is over the entire simulation since no transition occurs. The average reduction in <NHbonds>/molecule is maximum (0.1) when weak H-bonds are included (cutoff 3.3 Å, bottom row) and smallest (0.015) if only strong H-bonds (cutoff 2.8 Å, top row) are analyzed. In the 120 Å box (third column) and for the strong H-bonds the loss in <NHbonds>/molecule is insignificant but clearly increases when weak H-bonds are included. We also observe a significant decrease in the fluctuation of <NHbonds>/molecule between simulations in the smallest and the largest box sizes and a pronounced drop in <NHbonds>/molecule for the 75 and 90 Å boxes prior to the transitions, see insets.

https://doi.org/10.7554/eLife.35560.011

If the free energy difference between the T0 and R0states (7 kcal/mol, see above) were to arise entirely from the relative stability of the water-water network, this value corresponds to an energy difference of 104 kcal/mol per water molecule when comparing the 90 and 150 Å boxes, which differ by 80,000 water molecules. Such small energy differences are very difficult to capture in MD simulations and is not attempted here. However, it is interesting to note that the average number of hydrogen bonds per water molecule, <NHbonds>/molecule (see Figure 3) shows such an effect: for the three smaller boxes the <NHbonds>/molecule decreases by 103 to 104 with every transition. This is consistent with the estimated energy change per water molecule. Furthermore, Figure 3 demonstrates that the fluctuation of <NHbonds>/molecule decreases with increasing box size which is the behaviour expected from statistical mechanics. It should be pointed out that the running averages were evaluated over time intervals between which transitions were observed, see Figure 1.

Based on hydrodynamic arguments, Yeh and Hummer (2004) showed that the water self-diffusion coefficient, D, calculated from an MD simulation of pure water (without ions) in a periodic box, scales as N1/3, where N is the number of particles. For the largest box they studied (40 Å) the size correction was negligible. Interestingly, our simulation of the 75 Å box containing Hb yielded a value of D that is much too small (D=4.25×105 cm2/s vs D=5.95×105 cm2/s); the latter is the correct value for the TIP3P water model. In Figure 4A, we show the results for the value of D as a function of box size, plotted versus 1/L (nm1). As expected, all the pure water boxes are large enough so the calculated value agrees with the extrapolated value of Yeh and Hummer (D=5.9×105cm2/s) within statistical error. However, the calculated self-diffusion coefficient from MD simulations with Hb present, is identical to that of pure water only for the largest 150 Å box.

System-size dependence of the water diffusion coefficient.

(A) Water diffusion coefficients D as a function of system size for systems with and without protein. Also, results from Hummer et al. are shown for comparison and validation. (B) The water D as a function of time for the first instability to occur. For the 75, 90, and 120 Å boxes instabilities were observed (see Figure 1) and scale linearly with the water D. The yellow symbol for the 150 Å box is the extrapolated value (700 ns).

https://doi.org/10.7554/eLife.35560.013

The results described here suggest that the correct free energy of hemoglobin, at least to the extent that T0 is stable relative to R0, requires that the simulation be done in a box large enough so that the water environment behaves like bulk water. In Figure 4B, we plot the calculated Dvalue for boxes containing Hb versus the time in ns when the first transition from the T0 structure takes place. As can be seen, there is a linear relationship between the two. Of most interest is the fact that an extrapolation of the line indicates that in the 150 Å box, the first transition away from T0 should take place at 700 ns. However, we have shown that in the 150 Å box, T0 is still present at 1.4 μs. This provides strong evidence that in a 150 Å box, T0 is in fact stable. This linear dependence was not expected. It is an interesting result whose origin, though, still requires explanations at a molecular level. The result that the lifetime of the T0 state increases systematically with the increase in box size effectively corresponds to multiple simulations. The T-state was finally found to be stable in the 150 Å box for 1.4 μs, significantly longer than the extrapolated lifetime value (700 ns). The idea that μs-plus simulations are needed has become a ‘lore’ (similar to the box size-related ‘rule’ investigated here) with the availability of bigger computers, even when they are not required for a particular problem, as is the case here.

Simulations for all box sizes have also been done with His146 deprotonated. The results for that system showed that T0 is stabilized in larger boxes, but after less than 100 ns a transition to an R-like state occurs. This provides direct evidence that His146 protonation is essential for stabilizing T0, in accord with the Perutz model (Perutz, 1970).

To relate the above results to the hydrophobic effect, we use the construct of Chandler (Chandler, 2005), who showed (Figure 5A) that significant water depletion around a spherical hydrophobic solute is expected when its radius is larger than 1 nm (10 Å). Since Hb is nearly spherical with a 2.8 nm radius, we calculated the g(r) for the 90 and 150 Å box (see Figure 5B). The behavior in the 150 Å box is consistent with the expected water structuring, while for the 90 Å box it is not.

Figure 5 with 4 supplements see all
System-size dependence of the average equilibrium water density.

(A) From Ref. (Chandler, 2005) (Figure 3); the average equilibrium water density a distance r+R from spherical cavities. R is the (spherical) size of the solute. Solid lines for the ideal hydrophobe, dashed line when van der Waals interactions are present. (B) Results for Hb in water for the 90 Å (red) and 150 Å (green) box. The radius of Hb (2.8 nm) is intermediate between the two cases in panel A.

https://doi.org/10.7554/eLife.35560.014

The full radial distribution function is reported in Figure 5—figure supplement 1 and Figure 5—figure supplement 2B shows the number of water molecules derived from it. It is found that up to a distance of 8.5 Å from the central cylinder 150 water molecules are present which is consistent with explicit counting, see Figure 5—figure supplement 2A, and Figure 5—figure supplement 3 for the corresponding probability distribution function. Furthermore, the structural transitions are accompanied by pronounced dewetting and water penetration as shown in Figure 5—figure supplement 4 for the 90 Å box.

To summarize, the T0 state is only found to be thermodynamically stable if (i) the hydration water behaves like bulk water as judged from the self-diffusion coefficient, (ii) the number of hydrogen bonds per water molecule is large enough and its fluctuation around the average sufficiently small, see Figure 3—figure supplement 1. Therefore, if water is not engaged in an adequate number of water-water H-bonds, solvent water is prone to attack the protein salt bridges, destabilizing them and eventually to break them.

As a more local measure of the effect of undersolvation the radial distribution function g(r) around the (His146β1)CG–O(water) was determined for all box sizes, see Figure 6. For the 150 Å box the g(r) remains almost invariant whereas for the 120 Å box an appreciable change occurs during the time when the major structural transition at 920 ns takes place. For the smallest boxes, the local g(r) is very variable, which supports the importance of locally structured water molecules for stabilizing the T0 state.

Averaged radial distribution functions g(r) between His146 and water for the different box sizes and for simulation windows before, during and after the structural transitions (see Figure 1).

Analysis for (A) the C-terminal (CT) of His146 and Water H (Hw) and (B) for C of His146 and water O (Ow). For the largest (150 Å) box no appreciable change in g(r) is found. The two peaks at 2.7 and 4.1 Å in (A) and the two peaks at 5 and 7.5 Å in (B) are due to a water network as indicated in the sketches. The sketches in A and B represent different views of the same snapshot taken from the MD simulation of the 150 Å box at t=960 ns and is used here as reference to describe the water network because this box represented a stable g(r). Sketches in (A) emphasise the water network at the C-terminal of His146, and sketches in (B) emphasise the water network at the Cγ of His146. Waters W0, W1 and W2 are key structural waters (used as anchor) in the stability of the T state (they are framed with a black line). This network is stable for simulations in the largest box but unstable in the other boxes. His143 is singly protonated (Nϵ) and His146 is doubly protonated (Nδ and Nϵ).

https://doi.org/10.7554/eLife.35560.019

In addition to the global behavior, local structural changes involving the two His146 residues were examined. The results (Figure 1—figure supplement 35) show that for the largest box, the parameters examined (e.g. the salt bridge to Lys40, discussed by Perutz) fluctuate around the equilibrium values near those of the T0 structure while for the smaller boxes there are abrupt changes, correlated with global structural transitions. The transition away from the T0 structure extends over approximately 10 ns during which several important contacts are broken or formed, see Figure 1—figure supplement 5A. This value suggests that the apparent activation energy from the simulations in the smaller solvent boxes is near zero, in contrast to an estimate on the order of seconds for the transition if the barrier arose from the 7 kcal/mol difference in stability.

Importantly, the ability of the 150 Å box to stabilize the T-state starting from the R-like structure in the 120 Å box after 1 μs was also explored. Solvating this structure in the 150 Å box, minimizing, heating and equilibrating it (see Materials and methods) and following the equilibrium dynamics for 1 μs, the final RMSD compared to the 2DN2 (T-state) structure is 1.97 Å (2.55 compared to 2DN3 (R-state)) starting from an RMSD of 2.73 Å (with respect to 2DN2) and 1.42 Å (2DN3) at the beginning of the simulation (after heating and equilibration), respectively. Concomitantly, the His146β1–His146β2 and His143β1–His143β2 distances change from values typical for an R-state structure (both 12.6 Å) to 14.0 Å (for His143) and 18.0 (for His146) and approach separations indicative of a T-state structure (18.6 Å and 30.9 Å) without, however, fully completing the transition to a T-state structure. As mentioned above, the R0 to T0 transition time is about 20 μs, much longer than the simulation time.

Given the increasing use of molecular dynamics simulations to study conformational transitions in large proteins and in an explicit solvent environment, the present result that much larger boxes than those used standardly are required for what appears to be the hydrophobic effect to be manifested is of general interest. It has wide ranging implications for the interpretation and validity of previous simulations, as well as those to be undertaken in the future. Given that the magnitude of the effect is expected to depend on the size (and shape) of the molecule and its hydrophobicity, as well as possibly other properties (Chandler, 2005), the requirement for the use of larger boxes in simulations will have to be investigated in each case. One particularly relevant situation where capturing the correct diffusional dynamics of the environment will play a crucial role is that in atomistic simulations for the crowded conditions (Feig et al., 2017) that exist in entire cells or parts thereof (Zhou et al., 2008). Possible errors in molecular dynamics simulation results for other phenomena, such as protein folding, for example, as well as for polymeric materials more generally have to be considered as well.

Materials and methods

The influence of solvent layers on the structural stability of hemoglobin tetramer is investigated using Molecular Dynamics (MD) simulations. Extended large-scale simulation were performed with the CHARMM36 all atom force-field (Best et al., 2012) and the TIP3P water model (Jorgensen et al., 1983) using version 5.1.4 of the GROMACS package (Abraham et al., 2015) on GPUs.

The coordinates of the starting structure are taken from the X-ray structure of tetrameric human hemoglobin in the deoxy form, PDB code 2DN2 (1.25 Å resolution) (Park et al., 2006). The protonation states of the histidines were based on the 2013 study of Zheng et al. and the terminal β histidines (His146) were both doubly protonated (Zheng et al., 2013). The protein was solvated in four different cubic boxes of increasing size: 75, 90, 120 and 150 Å. The system was neutralized by adding counter ions and the salt concentration of 0.15 m/L was achieved using Na+ and Cl. The total number of atoms is: 39,432 for the 75 Å box, having 10,543 water molecules and 42 Na+ and 38 Cl ions; 72,142 for the 90 Å box, having 20,756 water molecules and 70 Na+ and 66 Cl ions; 163,480 for the 120 Å box, having 53,287 water molecules and 160 Na+ and 156 Cl ions; 318,911 for the 150 Å box, having 105,073 water molecules and 309 Na+ and 305 Cl ions.

For the electrostatic interactions, particle-mesh Ewald (PME) was used with a grid spacing of 1 Å, a relative tolerance of 106 and a cutoff of 10 Å, together with a 10 Å switching for the Lennard-Jones (LJ) interactions. The LINCS algorithm (Hess et al., 1997) was used for constraining bonds involving H-atoms. Each system was first energy minimized for 50,000 steps using steepest descent, heated from 0 to 300 K in increments of 10 K in NVT for 300 ps, followed by 500 ps (NVT) and 500 ps of NpT equilibration at p=1 atm with a time step of 2 fs. The center of mass of the protein was restrained to the center of mass of the simulation box. The Velocity Rescaling (Bussi et al., 2007) (with τ=0.1 ps) and Parrinello-Rahman (Parrinello and Rahman, 1981) methods were used for temperature and pressure control, respectively. The velocity rescaling method is an extension of the Berendsen thermostat to which a stochastic force chosen such as to generate a correct canonical distribution is added (Bussi et al., 2007). The MD simulations for all systems were carried out for at least 1 μs at constant temperature and pressure (NpT) at 300 K and one atm with a time step of 2 fs and the 1σ temperature fluctuations over the 1 μs trajectories were 0.1 K for the 150 Å box and 0.5 K for the 90 Å box.

Water self-diffusion coefficients D were calculated for box sizes 75, 90, 120 and 150 Å, in the presence and in the absence of Hb, over the entire 1 μs trajectory. In the absence of the protein, the simulation boxes contained pure water systems (no ions were included). Including ions at physiological concentrations will typically change the water self-diffusion coefficient by 1% to 2% (Kim et al., 2012). First, the mean square displacement (MSD) of all oxygen atoms from a set of initial positions was calculated using mass-weighted averages. Then, the diffusion constant was calculated from the slope of the mean-squared displacement DPBC=Limtt|r(t)r(0)|26 averaged over all water molecules of a particular trajectory. Errors are estimated from the difference of the diffusion coefficients obtained from separate fits over the two halves of the fit interval.

The TIP3P self diffusion coefficient calculated in the simulation by Yeh et al. (using periodic boundary conditions) is DPBC=5.8×105 cm2/s (Yeh and Hummer, 2004). However, simulations of water at ambient conditions and a Lennard-Jones (LJ) fluid show that the diffusion coefficients depend on system size (Allen and Tildesley, 1987). Thus, the diffusion coefficient corrected for system-size effects is D0=6.1×105 cm2/s for an infinite system of TIP3P water at 298 K and ambient temperature (Yeh and Hummer, 2004). For direct comparison of our values the error κ=L×(D0DPBC) was subtracted from D0. The resulting TIP3P value is D0,corr=5.95×105 cm2/s.

In order to directly compare with the literature (Yeh and Hummer, 2004) the 40 Å box was also considered here. The literature value of DPBC=5.65±0.16 (105 cm2/s) compares with DPBC=5.59±0.013 (105 cm2/s) computed here which validates the present simulations.

Appendix

Structural analyses

The distance between the two Cα atoms of His146 of the two β chains is one meaningful coordinate to distinguish the T and the R configuration. The distance between the two Cα atoms of His143 present at the opening of the central cavity can be used to monitor the width of the opening. The RMSD between the Cαs relative to the 2DN2-Xray structure was also evaluated.

For the 75 Å box the collapse after 130 ns leads to a concomitant increase in the Cα-Cα RMSD from 2DN2 by more than 1 Å, a repositioning of the β2 chain relative to the other three chains and a decrease in the degree of hydration of the central channel by up to 20%. Similar observations apply to the 90 (120) Å boxes for which three structural transitions are observed: one after 480 (630) ns (Step 1), the next one after 780 (820) ns (Step 2), and Step 3 after 880 (920) ns. The transitions in the smaller boxes involve different changes (see Figure 1—figure supplement 2). Also, the end point structures after 1 μs are not the same (see Table 1).

In order to evaluate and monitor the opening of the central cavity the angle between α1β1 and α2β2 sub-units groups is calculated as the angle formed between their planes (Figure 1—figure supplement 2). The plane is defined by three atoms: Cα of His146β, Feβ and Feα. The smaller the angle the closer are the subunits and the tighter is the central cavity. This affects the number of water molecules that can access the channel. The water count was thus evaluated and confirmed this observation.

The number of water molecules inside the central cylinder of the protein was also analyzed see Figure 5—figure supplements 3 and 4. The water count includes only waters in the central cylinder, which was defined as a cylinder of 8.5 Å base area, 22 Å of height and having the four heme moieties as the center of mass, and not waters at the tetramerization interface.

The ion distribution was considered, as well. The concentration of ions was 0.15 mol/L for all box sizes. It is known that the ionic strength of the solution plays a role in the oxygenation behaviour. (Szabo and Karplus, 1972) Also, the local ion concentration has an indirect influence by modulating the water activity, which in turn affects allosteric regulation of Hb. (Salvay et al., 2003) The ion diffusivity changes from 2.1 to 2.8×105 cm2/s for simulations in the 75 to 150 Å boxes which also reflects the increased diffusivity of the bulk water molecules with increasing box size, see Figure 4A. The average number of ions within a 20 Å sphere around the center of mass of Hb is 9.6 (75 Å), 8.5 (90 Å), 8.7 (120 Å), and 9.8 (150 Å). Hence, the ion distribution does not drive the transition directly, which is consistent with additional test simulations carried out, as described below.

Additional modifications considered

Further modifications for which the stability of the T-state was investigated included (i) strategic placements (by using constraints) of 4 Cl ions at the β1β2 interface where the ligand 2,3-diphosphoglycerate (DPG) binds, (ii) reduced force constants for the Fe-NHis stretching to render the heme group more adaptable to global structural changes of the protein, (iii) placement of the DPG ligand which is known to bind to deoxy-Hb, (iv) modified charges of the hydrogen bonds to stabilize the hinge contacts (Jones et al., 2014), (v) different protonation state of His146β1 and His146β2 (δ, ϵ, and doubly protonated) and (vi) scaling the water-protein nonbonded interactions, specifically the van der Waals well depths. (Best et al., 2014) All these simulations were run starting from the T-state in the 75 Å box and were found to become unstable within 70 to 200 ns. The outcome of the last two modifications considered (v and vi) is reported in Figure 1—figure supplement 6. The figure shows the Cα–Cα distance between residues His146β1 and His146β2 for the T-state Hb in the 75 Å box during 400 ns of 4 MD simulations for the variants (1) protonated His146 (protonated Nδ and Nϵ) with (black) and (2) without (red) scaled protein-water interactions; (3) unprotonated His146 (proton on Nϵ green) and (4) unprotonated His146 (proton on Nδ blue) with scaled protein-water interactions. In all cases the T-state structure becomes unstable to various degrees and decays more or less completely towards the R-state. The red trace in Figure 1—figure supplement 6 is identical to the black trace in Figure 1 of the main manuscript.

Particular observations made for these additional simulations are as follows: In the presence 4 Cl ions in the cavity where DPG binds the T-state is stable for 45 ns. Constraining the chloride ions does not significantly extend the time during which the protein is in its T-state. The rationale behind reducing the force constants for the Fe-N bonds is that by softening the internal degrees of freedom the protein as a whole can adapt more easily to perturbations induced by strain. However, even reducing the stretching force constant by a factor of 10 only increases the lifetime of the T-state out to 60 ns. The DPG ligand was placed at the β1β2 interface and harmonic constraints were applied to keep it close to this initial position. Despite the present of the ligand, the Cα–Cα distance between His146β1 and His146β2 decreases to 20 Å within 20 ns and hence the T-state decays. Recent experimental work (Jones et al., 2014) suggested that special hydrogen bond contacts (so called ‘hinge’ contacts) exist at the α1β2 interface and their presence and absence have been related to the TR transition. In order to stabilize these contacts (Tyr42–Asp99 and Trp37–Asp94), the partial charges of the hydrogen and acceptor atoms were increased by 10 and 20%, respectively. With this, the T-state lifetime was 50 ns or less, compared to 40 ns if the charges were unmodified. Such tests had also been done for the T-state by constraining the Perutz H-bonds for 100 ns in equilibrium simulations. (Hub et al., 2010) However, after removing the constraints, the T-structure decayed again during the following 100 ns simulation. Finally, it had been found that a strengthening of the protein-water van der Waals interaction by 10% is sufficient to recover the correct dimensions of intrinsically disordered and unfolded proteins. (Best et al., 2014) Applying this to the apparent instability of the T-state of Hb does, however, not extend its lifetime. After 75 to 120 ns the T-state decays towards the R-state structure (see Figure 1—figure supplement 6). We also note that, although in the initial simulation of the 75 Å box with His 146 doubly protonated, the Cα-Cα distance between the His 146 residues decreases by only a small amount (from 31 Å to 25 Å in step 1), by other criteria the decaying T-state approaches the R-state (see Figure 1—figure supplement 2).

References

  1. 1
  2. 2
    Computer Simulation of Liquids
    1. MP Allen
    2. DJ Tildesley
    (1987)
    Oxford, U.K: Clarendon Press.
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
    Proceedings of the International School of Physics Enrico Fermi
    1. D Chandler
    2. P Varilly
    (2012)
    75–111, Lectures on molecular- and nano-scale fluctuations in water, Proceedings of the International School of Physics Enrico Fermi, 176.
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
    Quaternary conformational changes in human hemoglobin studied by laser photolysis of carboxyhemoglobin
    1. CA Sawicki
    2. QH Gibson
    (1976)
    The Journal of Biological Chemistry 251:1533–1542.
  31. 31
  32. 32
  33. 33
    Observation of the dissociation of unliganded hemoglobin
    1. JO Thomas
    2. SJ Edelstein
    (1972)
    The Journal of Biological Chemistry 247:7870–7874.
  34. 34
  35. 35
  36. 36
  37. 37

Decision letter

  1. Yibing Shan
    Reviewing Editor; DE Shaw Research, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Valid Molecular Dynamics Simulations of Human Hemoglobin Require a Surprisingly Large Box Size" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by John Kuriyan as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

This is a potentially interesting study in which the authors suggest that the stability of the T state of hemoglobin in molecular dynamics (MD) simulations depends critically on the size of the simulation box and that the T state remains stable only with a sufficiently large box containing a relatively large portion of bulk-like solvent. While this observation may bear important practical implications to MD simulation of biomolecules (and polymers in general), its statistical validity needs to be strengthened with further data, and the understanding of its origin requires further analyses.

Essential revisions:

1) The key observation was made based on four simulations in four different simulation boxes. Each simulation has been performed only once. Given the stochastic nature of protein conformational change in MD simulations, repeats of these simulations are essential to show that the observation is of statistical significance. Minimally, each simulation should be repeated once or twice, which should be practically manageable since each simulation is 1 microsecond long. Other simulations should be considered to help establish the key observation. For example, in the cover letter, the authors stated that they "ran simulations starting with a partially decayed (i.e. R-state-like) structure from a 120 ˚A3 box in the 150 ˚A3 box and demonstrated that it progressed toward the T0 structure." These data would be very important, but it is not shown in the manuscript. The authors should show plots of all the distance, angles and RMSDs for this simulation so that it would be easier to assess the significance of their finding. They should also run the same system for an additional 1μs in the 120A box, to provide a fair comparison between the two. Various attempts by the authors to stabilize the T state, which is mentioned in the cover letter but largely not discussed in the main text or the supplemental information, should be reported with some details as well.

2) A possible reason for the instability of the T state would be inaccuracies in the potential energy function, in particular for the Heme group. Also, the 3-point models used to describe water in the simulations are known to provide only a rather rough description of liquid water properties and there is no reason to believe that they would accurately capture the subtle balance of electrostatic and hydrophobic forces that modulate the relative stability of the T and R states. Comparison with previous works is complicated by the use of a different force-field (CHARMM36 (?), instead of the Amber and the Gromos force-fields used in previous works). These issues need to be discussed. Also, the reference for the Charmm36 force-field, which refers to paper for parametrization of the Charmm lipid force-field by mistake, needs to be corrected.

3) The analyses to understand the origin of the possible phenomenon consisted of three parts, respectively concerning hydrogen bonds, diffusion, and radial distribution function of water molecules. The first two analyses showed that, with a larger simulation box, an average water molecule in the simulation increasingly takes up the bulk-water behavior. This is certainly sensible but also to be expected, as in a large box a larger portion of waters are not in the vicinity of the protein. Hemoglobin is a protein that contains a number of cavities and a water-filled channel. It is expected that water molecules in the cavities or in the channel will have a smaller diffusion coefficient and a smaller number of hydrogen bonds with respect to the bulk water. The analysis of the average water diffusion coefficient and number of hydrogen bonds as a function of box size (Figure 3 and 4) is merely showing convergence towards the bulk as the ratio between (water in cavities)/(bulk water) becomes smaller and smaller. The authors need to more clearly state their connection with the T stability of hemoglobin.

4) The authors' underlying hypothesis appears to be that the deviation of water molecules in small boxes from the bulk behavior affects the magnitude of the hydrophobic effect, and that this effect differentially impacts the stability of the two states because the T state is structurally more compact and with less solvent exposed area than the R state. If so, this reasoning should be stated explicitly with necessary data.

5) Figure 5A shows that the water RDF around a spherical solute approaches bulk behavior ~0.5-1 nm away from the solute, suggesting that a periodic box that is ~2 nm larger than the solute should give acceptable result. Figure 5B, however, is apparently at odds with the g(r) presented in Figure 5A, as bulk water behavior seems to be achieved only ~2.5 nm away from Hemoglobin. Figure 5—figure supplement 1 presents the whole RDF profile. This RDF calculated in this manuscript appears to be wrong, as there is water in the center of Hemoglobin, so the RDF should have value close to 1 in the center. If this is a mistake, it should be corrected. Did the authors checked that at the beginning of the simulation water was actually placed in the central cavities of the protein and that this water was properly equilibrated?

https://doi.org/10.7554/eLife.35560.023

Author response

Essential revisions:

1) The key observation was made based on four simulations in four different simulation boxes. Each simulation has been performed only once. Given the stochastic nature of protein conformational change in MD simulations, repeats of these simulations are essential to show that the observation is of statistical significance. Minimally, each simulation should be repeated once or twice, which should be practically manageable since each simulation is 1 microsecond long. Other simulations should be considered to help establish the key observation. For example, in the cover letter, the authors stated that they "ran simulations starting with a partially decayed (i.e. R-state-like) structure from a 120 ˚A3 box in the 150 ˚A3 box and demonstrated that it progressed toward the T0 structure." These data would be very important, but it is not shown in the manuscript. The authors should show plots of all the distance, angles and RMSDs for this simulation so that it would be easier to assess the significance of their finding. They should also run the same system for an additional 1μs in the 120A box, to provide a fair comparison between the two. Various attempts by the authors to stabilize the T state, which is mentioned in the cover letter but largely not discussed in the main text or the supplemental information, should be reported with some details as well.

As agreed with the Editor we did not repeat these simulations.

2) A possible reason for the instability of the T state would be inaccuracies in the potential energy function, in particular for the Heme group. Also, the 3-point models used to describe water in the simulations are known to provide only a rather rough description of liquid water properties and there is no reason to believe that they would accurately capture the subtle balance of electrostatic and hydrophobic forces that modulate the relative stability of the T and R states. Comparison with previous works is complicated by the use of a different force-field (CHARMM36 (?), instead of the Amber and the Gromos force-fields used in previous works). These issues need to be discussed. Also, the reference for the Charmm36 force-field, which refers to paper for parametrization of the Charmm lipid force-field by mistake, needs to be corrected.

Given that other simulations (Hub et al., 2010 and Yussuf et al., 2012) using different MD codes and force fields find similar instabilities it is unlikely that the FF parametrizations are the prime reason for the thermodynamic instability, see main text, third paragraph. Also, some of the additional tests that were carried out included FF-modifications which, however, did not lead to stabilization in the small water boxes. The reference to the CHARMM36 force field was corrected.

3) The analyses to understand the origin of the possible phenomenon consisted of three parts, respectively concerning hydrogen bonds, diffusion, and radial distribution function of water molecules. The first two analyses showed that, with a larger simulation box, an average water molecule in the simulation increasingly takes up the bulk-water behavior. This is certainly sensible but also to be expected, as in a large box a larger portion of waters are not in the vicinity of the protein. Hemoglobin is a protein that contains a number of cavities and a water-filled channel. It is expected that water molecules in the cavities or in the channel will have a smaller diffusion coefficient and a smaller number of hydrogen bonds with respect to the bulk water. The analysis of the average water diffusion coefficient and number of hydrogen bonds as a function of box size (Figure 3 and 4) is merely showing convergence towards the bulk as the ratio between (water in cavities)/(bulk water) becomes smaller and smaller. The authors need to more clearly state their connection with the T stability of hemoglobin.

In the sixth paragraph of the Results and Discussion it is now clarified that if the water self diffusivity and the hydration network resemble that of bulk water, T0 is stable. Figure 3—figure supplement 1 was added to substantiate this. The shape of the radial distribution functions is rather a consequence than the origin of the T0 stability.

4) The authors' underlying hypothesis appears to be that the deviation of water molecules in small boxes from the bulk behavior affects the magnitude of the hydrophobic effect, and that this effect differentially impacts the stability of the two states because the T state is structurally more compact and with less solvent exposed area than the R state. If so, this reasoning should be stated explicitly with necessary data.

We have added a statement to the first paragraph of the Results and Discussion which makes clear what is in the papers of Chothia et al. (1976) and Lesk et al. (1985) concerning the role of the hydrophobic effect in stabilizing the T0 state. What we do observe is, that the T0 state is thermodynamically stable if the water box is sufficiently large such that solvent water behaves as bulk water as judged from the self diffusion coefficient (see Figure 4) and the hydrogen bonded network (see new Figure 3—figure supplement 1).

5) Figure 5A shows that the water RDF around a spherical solute approaches bulk behavior ~0.5-1 nm away from the solute, suggesting that a periodic box that is ~2 nm larger than the solute should give acceptable result. Figure 5B, however, is apparently at odds with the g(r) presented in Figure 5A, as bulk water behavior seems to be achieved only ~2.5 nm away from Hemoglobin. Figure 5—figure supplement 1 presents the whole RDF profile. This RDF calculated in this manuscript appears to be wrong, as there is water in the center of Hemoglobin, so the RDF should have value close to 1 in the center. If this is a mistake, it should be corrected. Did the authors checked that at the beginning of the simulation water was actually placed in the central cavities of the protein and that this water was properly equilibrated?

There are ∼ 130 water molecules in the central cavity of Hb, see new Figure 5—figure supplements 2 and 3. This is from two different analyses. One involves explicit counting and the other one integration of g(r). The time series confirm that the central cavity is filled and that the two counts are consistent.

https://doi.org/10.7554/eLife.35560.024

Article and author information

Author details

  1. Krystel El Hage

    Department of Chemistry, University of Basel, Basel, Switzerland
    Contribution
    Conceptualization, Data curation, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4837-3888
  2. Florent Hédin

    Department of Chemistry, University of Basel, Basel, Switzerland
    Contribution
    Data curation, Formal analysis, Investigation, Writing—original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6341-7557
  3. Prashant K Gupta

    Department of Chemistry, University of Basel, Basel, Switzerland
    Contribution
    Data curation, Formal analysis, Investigation, Writing—original draft
    Competing interests
    No competing interests declared
  4. Markus Meuwly

    Department of Chemistry, University of Basel, Basel, Switzerland
    Contribution
    Conceptualization, Resources, Software, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    m.meuwly@unibas.ch
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7930-8806
  5. Martin Karplus

    1. Department of Chemistry and Chemical Biology, Harvard University, Cambridge, Massachusetts, United States
    2. Laboratoire de Chimie Biophysique, ISIS, Université de Strasbourg, Strasbourg, France
    Contribution
    Conceptualization, Investigation, Methodology, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    marci@tammy.harvard.edu
    Competing interests
    No competing interests declared

Funding

Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (NCCR MUST)

  • Markus Meuwly

Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (200020-169079)

  • Markus Meuwly

Charmm Development Project

  • Martin Karplus

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank David Chandler, to whom we dedicate this paper, for many fruitful discussions. The work in Switzerland was supported by the Swiss National Science Foundation through grants 200021–117810, and the NCCR MUST. The work at Harvard was supported in part by the CHARMM Development Project.

Reviewing Editor

  1. Yibing Shan, DE Shaw Research, United States

Publication history

  1. Received: February 1, 2018
  2. Accepted: June 18, 2018
  3. Version of Record published: July 12, 2018 (version 1)

Copyright

© 2018, El Hage et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,137
    Page views
  • 145
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Structural Biology and Molecular Biophysics
    Kouki K Touhara, Roderick MacKinnon
    Research Article
    1. Structural Biology and Molecular Biophysics
    Fei Ye et al.
    Research Article