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

Recent molecular dynamics (MD) simulations of human hemoglobin (Hb) give results in disagreement with experiment. Although it is known that the unliganded (T[Formula: see text]) and liganded (R[Formula: see text]) tetramers are stable in solution, the published MD simulations of T[Formula: see text] undergo a rapid quaternary transition to an R-like structure. We show that T[Formula: see text] 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 T[Formula: see text] tetramer, to be manifested. Even in the largest box, T[Formula: see text] 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.

Introduction

Human hemoglobin is the paradigmatic model system for n class="Chemical">cooperativity in proteins. It transports oxygen from the lungs to the tissues and it is composed of two identical -chains ( and ) and two identical chains ( and ). They form two identical dimers ( and ), whose relative orientation differs significantly in the unliganded (T) and liganded (R) 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 (T) and oxy (R)) 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 T and R for which 1.25 Å resolution X-ray structures are available (T, 2DN2; R, 2DN3) (Park et al., 2006). Although R was found to be stable for several hundred nanoseconds, the T 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 T state is about seven kcal/mol more stable than the R state (Edelstein, 1971), which is derived from the ratio of the dissociation constants of liganded and unliganded Hb of  (Thomas and Edelstein, 1972).

Figure 1—figure supplement 2.

Global structural changes depending on box size.

Panel A: Overall structure of the hemoglobin tetramer. The iron atoms (Fe, yellow spheres) and His146 side chains (green spheres) are specifically indicated. The angle measures the angle between the two planes containing His146FeFe and His146FeFe, respectively. Panel B: temporal change of (from top to bottom) the C–C distance between His146 and His146, the C–C distance between His143 and His143, C RMSD relative to the 2DN2 X-ray structure, and the angle . The dashed lines (black, orange, red) indicate transition points for the 75, 90, and 120 Å boxes, respectively. Cyan and blue arrows indicate the values of the corresponding observables found for the deoxy T (2DN2) and oxy R (2DN3) states, respectively. Top panel is also given in the main MS.

Although the rate of the T to R transition has not been measured directly, it can be estimated from the experimentally determined R to T transition rates and the R/T equilibrium constant. For the unliganded hemoglobin tetramer, the R to T rate is about 20 s (Sawicki and Gibson, 1976). With the equilibrium constant of [T/R] equal to  (Edelstein, 1971; Thomas and Edelstein, 1972), the T to R 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 T and R. 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 T 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 T 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).

Tetrameric hemoglobin solvated in boxes of size: (from left to right) 75, 90, 120 and 150 Å. Spheres correspond to Na and Cl ions, added to neutralize the system and to achieve a biologically appropriate n class="Chemical">salt concentration of 0.15 mol/L.

From top to bottom: the Ramachandran -angle of His143, n class="Chemical">Lys144, Tyr145 and His146. Left panels report time series for , and right panels those for . Cyan and blue arrows indicate the corresponding values found in crystal structure of the deoxy (2DN2) and oxy (2DN3) state, respectively. Color code as in Figure 1—figure supplement 1.

Figure 1—figure supplement 1.

Tetrameric hemoglobin solvated in boxes of different sizes.

Tetrameric hemoglobin solvated in boxes of size: (from left to right) 75, 90, 120 and 150 Å. Spheres correspond to Na and Cl ions, added to neutralize the system and to achieve a biologically appropriate salt concentration of 0.15 mol/L.

Local structural changes around His146.

Temporal structural changes of the C–C distance between His146 and His146 in the 75 Å box for different simulation conditions.

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 n class="Chemical">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 T tetramer dimensions are approximately Å, 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 ensemble (see SI for details).

To investigate the possibility that larger boxes were required for stabilizing T, we carried out 1s 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 and His146, 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 T 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 R structure (2DN3) for the 90 and 120 Å simulations and the deoxy T structure (2DN2) for the 150 Å box simulation (see also Table 1).

Figure 1.

Global structural changes depending on box size.

Panel A: Overall structure of the hemoglobin tetramer. The His146 side chains (green spheres) are specifically indicated. Panel B: Temporal change of the C–C distance between His146 and His146. 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 T (2DN2) and oxy R (2DN3) states, respectively.

Tetrameric hemoglobin solvated in boxes of size: (from left to right) 75, 90, 120 and 150 Å. Spheres correspond to Na and Cl ions, added to neutralize the system and to achieve a biologically appropriate salt concentration of 0.15 mol/L.

Panel A: Overall structure of the hemoglobin tetramer. The iron atoms (Fe, yellow spheres) and His146 side chains (green spheres) are specifically indicated. The angle measures the angle between the two planes containing His146FeFe and His146FeFe, respectively. Panel B: temporal change of (from top to bottom) the C–C distance between His146 and His146, the C–C distance between His143 and His143, C RMSD relative to the 2DN2 X-ray structure, and the angle . The dashed lines (black, orange, red) indicate transition points for the 75, 90, and 120 Å boxes, respectively. Cyan and blue arrows indicate the values of the corresponding observables found for the deoxy T (2DN2) and oxy R (2DN3) states, respectively. Top panel is also given in the main MS.

From top to bottom: the Ramachandran -angle of His143, Lys144, Tyr145 and His146. Left panels report time series for , and right panels those for . Cyan and blue arrows indicate the corresponding values found in crystal structure of the deoxy (2DN2) and oxy (2DN3) state, respectively. Color code as in Figure 1—figure supplement 1.

Left: Interactions involving His146. Right: Interactions involving His146. From top to bottom: [a] water-mediated contact between (His146)COO–OC(Pro37), [b] the salt bridge between (His146)COO–NZ(Lys40), [c] the contact between (His146)COO–NZ(Lys132), [d] the contact between (His146)COO and NE(His2) and [e] the salt bridge between (His146)NE2 and COO(Asp94). Cyan and blue arrows indicate the corresponding values found in crystal structure of the deoxy (2DN2) and oxy (2DN3) state, respectively. The dashed lines indicate breaking points along the 1s MD simulation and the black arrows indicates the first to break. The straight green line represents the averaged distance over 1 s of MD simulation for the 150 Å box as reference.

(A) Close-up of the [600;650] and [440;520] ns interval of Figure 1—figure supplement 1, illustrating the order of events for transitions in the 120 (left panel) and 90 Å box (right panel), respectively. Black arrows indicates the transmission of the structural modifications. (B) Structural details for every quantity from the 90 Å box at 10 and 800 ns, respectively. The dashed black arrow indicates the separation distance between and .

Temporal structural changes of the C–C distance between His146 and His146, depending on the protonation state of His 146 in the 75 Å box with and without scaled protein-water interactions.

Figure 2.

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 ns. (B) Superposition of the 2DN3 structure and the Hb structure in the 120 Å box at ns. (C) Superposition of the 2DN2 structure and the structure in the 150 Å box at ns. In each case the RMSDs (in Å) to the 2DN2 and 2DN3 are also given.

Table 1.

C RMSD (in Å) relative to the 2DN2 (T) and 2DN3 (R) 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.

Structure	C_α-Cα RMSD to 2DN2	C_α-Cα RMSD to 2DN3
2DN2	−	2.43
2DN3	2.43	−
75 Å box	2.37	2.59
90 Å box	3.35	0.64
120 Å box	2.44	0.30
150 Å box	0.39	3.09

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

If the free energy difference between the T and Rstates (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 kcal/mol per water molecule when comparing the 90 and 150 Å boxes, which differ by 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, /molecule (see Figure 3) shows such an effect: for the three smaller boxes the /molecule decreases by to with every transition. This is consistent with the estimated energy change per water molecule. Furthermore, Figure 3 demonstrates that the fluctuation of /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.

Figure 3.

The average number of H-bonds per water molecule /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 /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 /molecule is insignificant but clearly increases when weak H-bonds are included. We also observe a significant decrease in the fluctuation of /molecule between simulations in the smallest and the largest box sizes and a pronounced drop in /molecule for the 75 and 90 Å boxes prior to the transitions, see insets.

Left column: Probability distribution of number of H-bonds relative to the average of the 150 Å box. Right column: probability distribution of the fluctuation for each box individually. The left column illustrates the loss in the average number of relative to the 150 Å box for smaller box sizes. This loss increases with including medium ( Å) and weak ( Å) H-bonds. Solid lines for simulations including the protein, broken lines for pure water boxes. As can be seen (left column), for a box size of 150 Å (solid and broken green lines) the average number of Hbonds with and without protein changes only marginally, whereas for a 75 Å box these numbers (solid and broken black lines) differ considerably. Also, the depends sensitively on the box size.

Based on hydrodynamic arguments, Yeh and Hummer (2004) showed that the water self-diffusion coefficient, , calculated from an MD simulation of pure water (without ions) in a periodic box, scales as , where 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 that is much too small ( cm/s vs cm/s); the latter is the correct value for the TIP3P water model. In Figure 4A, we show the results for the value of as a function of box size, plotted versus (nm). As expected, all the pure water boxes are large enough so the calculated value agrees with the extrapolated value of Yeh and Hummer (cm/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.

Figure 4.

System-size dependence of the water diffusion coefficient.

(A) Water diffusion coefficients 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 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 . The yellow symbol for the 150 Å box is the extrapolated value (700 ns).

(A) Water diffusion coefficients 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 n class="Chemical">water 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 . The yellow symbol for the 150 Å box is the extrapolated value (700 ns).

The results described here suggest that the correct free energy of hemoglobin, at least to the extent that T is stable relative to R, requires that the simulation be done in a box large enough so that the water envn class="Chemical">ironment behaves like bulk water. In Figure 4B, we plot the calculated value for boxes containing Hb versus the time in ns when the first transition from the T 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 T should take place at 700 ns. However, we have shown that in the 150 Å box, T is still present at 1.4 s. This provides strong evidence that in a 150 Å box, T 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 T 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 T is stabilized in larger boxes, but after less than 100 ns a transition to an R-like state occurs. Tn class="Chemical">his provides direct evidence that His146 protonation is essential for stabilizing T, 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 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.

System-size dependence of the average equilibrium water density.

() From Ref. (Chandler, 2005) (Figure 3); the average equilibrium water density a distance from spherical cavities. 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.

A: From Ref.(Chandler, 2005) (Figure 3); the average equilibrium water density a distance from spherical cavities. is the (spherical) size of the solute and is the distance from the surface of the hydrophobe. Solid lines for the ideal hydrophobe, dashed line when van der Waals interactions are present. The radial distribution function reaches the bulk value Å away from the solute’s surface. B: results for Hb in water for the 90 Å (red) and 150 Å (green) box. The radius of Hb (2.8 nm or 28 Å) is intermediate between the two cases in panel A. The coordinate is split in two segments the first (running from 0 to 28 Å) is the region inside Hb. The second (from 0 to 45 Å) is the region outside Hb. The radial distribution function reaches the bulk value Å away from the solute’s surface.

Panel A: Time series for the number of water molecules in the central cylinder (see also Figure 5—figure supplement 3) for 90 Å and 150 Å boxes. Panel B: radial distribution function (solid lines) and (dashed lines) for 90 and 150 Å boxes. The left panel shows the entire range of distances from the center of the protein and the right panel only that concerning the central cylinder. The value of is 150 which is consistent with the results from explicit counting and confirms that the is meaningful.

Probability distribution of the number of water molecules inside Hb over 1 s of MD simulation depending on box size. The water count includes only waters in the central cylinder and not waters at the tetramerization interface. The region was defined as a cylinder of 8.5 Å radius and 22 Å of height and having as center the center of mass of the 4 Heme moieties.

The number of water molecules at the interface from 1 s MD for the 90 Å box. (A) Time series of the interface. (B) Time series of the interface. (C) Illustration of the studied regions. The transitions occur at 480, 780 and 880 ns, respectively.

() From Ref. (Chandler, 2005) (Figure 3); the average equilibrium water density a distance from spherical cavities. 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 n class="Chemical">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.

Figure 5—figure supplement 3.

Distribution of water molecules in the central cylinder for the different box sizes.

Probability distribution of the number of water molecules inside Hb over 1 s of MD simulation depending on box size. The water count includes only waters in the central cylinder and not waters at the tetramerization interface. The region was defined as a cylinder of 8.5 Å radius and 22 Å of height and having as center the center of mass of the 4 Heme moieties.

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 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.

Figure 5—figure supplement 1.

Enhanced version of Figure 5 from the main MS.

A: From Ref.(Chandler, 2005) (Figure 3); the average equilibrium water density a distance from spherical cavities. is the (spherical) size of the solute and is the distance from the surface of the hydrophobe. Solid lines for the ideal hydrophobe, dashed line when van der Waals interactions are present. The radial distribution function reaches the bulk value Å away from the solute’s surface. B: results for Hb in water for the 90 Å (red) and 150 Å (green) box. The radius of Hb (2.8 nm or 28 Å) is intermediate between the two cases in panel A. The coordinate is split in two segments the first (running from 0 to 28 Å) is the region inside Hb. The second (from 0 to 45 Å) is the region outside Hb. The radial distribution function reaches the bulk value Å away from the solute’s surface.

Figure 5—figure supplement 2.

Number of water molecules in the central cylinder of Hb.

Panel A: Time series for the number of water molecules in the central cylinder (see also Figure 5—figure supplement 3) for 90 Å and 150 Å boxes. Panel B: radial distribution function (solid lines) and (dashed lines) for 90 and 150 Å boxes. The left panel shows the entire range of distances from the center of the protein and the right panel only that concerning the central cylinder. The value of is 150 which is consistent with the results from explicit counting and confirms that the is meaningful.

Figure 5—figure supplement 4.

Water at the interface.

The number of water molecules at the interface from 1 s MD for the 90 Å box. (A) Time series of the interface. (B) Time series of the interface. (C) Illustration of the studied regions. The transitions occur at 480, 780 and 880 ns, respectively.

To summarize, the T 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.

Figure 3—figure supplement 1.

Number of H-bonds and their fluctuations for pure water and simulations including the protein for the four box sizes.

Left column: Probability distribution of number of H-bonds relative to the average of the 150 Å box. Right column: probability distribution of the fluctuation for each box individually. The left column illustrates the loss in the average number of relative to the 150 Å box for smaller box sizes. This loss increases with including medium ( Å) and weak ( Å) H-bonds. Solid lines for simulations including the protein, broken lines for pure water boxes. As can be seen (left column), for a box size of 150 Å (solid and broken green lines) the average number of Hbonds with and without protein changes only marginally, whereas for a 75 Å box these numbers (solid and broken black lines) differ considerably. Also, the depends sensitively on the box size.

As a more local measure of the effect of undersolvation the radial distribution function around the (His146)CG–O(water) was determined for all box sizes, see Figure 6. For the 150 Å box the 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 is very variable, which supports the importance of locally structured water molecules for stabilizing the T state.

Figure 6.

Averaged radial distribution functions 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 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 ns and is used here as reference to describe the water network because this box represented a stable . 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 () and His146 is doubly protonated (N and N).

In addition to the global behavior, local structural changes involving the two His146 residues were examined. The results (Figure 1—figure supplement 3–5) 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 T structure while for the smaller boxes there are abrupt changes, correlated with global structural transitions. The transition away from the T 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.

Figure 1—figure supplement 3.

Local structural rearrangement of the C-terminus of the chain.

From top to bottom: the Ramachandran -angle of His143, Lys144, Tyr145 and His146. Left panels report time series for , and right panels those for . Cyan and blue arrows indicate the corresponding values found in crystal structure of the deoxy (2DN2) and oxy (2DN3) state, respectively. Color code as in Figure 1—figure supplement 1.

Figure 1—figure supplement 5.

Structural transition details.

(A) Close-up of the [600;650] and [440;520] ns interval of Figure 1—figure supplement 1, illustrating the order of events for transitions in the 120 (left panel) and 90 Å box (right panel), respectively. Black arrows indicates the transmission of the structural modifications. (B) Structural details for every quantity from the 90 Å box at 10 and 800 ns, respectively. The dashed black arrow indicates the separation distance between and .

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–His146 and His143–His143 distances change from values typical for an R-state structure (both Å) 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 R to T 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 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 for 300 ps, followed by 500 ps () and 500 ps of equilibration at 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 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 at 300 K and one atm with a time step of 2 fs and the 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 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 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 cm/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 efn class="Chemical">fects is cm/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 was subtracted from . The resulting TIP3P value is cm/s.

In order to directly compare with the literature (Yeh and Hummer, 2004) the 40 Å box was also considered here. The literature value of ( cm/s) compares with ( cm/s) computed here which validates the present simulations.

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 Rn class="Disease">MSDs 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 afn class="Chemical">fects 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?

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 Rn class="Disease">MSDs 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 n class="Chemical">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 n class="Chemical">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 afn class="Chemical">fects 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.