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
Cite this article as: eLife 2018;7:e35560 doi: 10.7554/eLife.35560
6 figures, 1 table and 1 additional file

Figures

Figure 1 with 6 supplements
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
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.

https://doi.org/10.7554/eLife.35560.003
Figure 1—figure supplement 2
Global structural changes depending on box size.

Panel A: Overall structure of the α1β1α2β2 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 His146β1Feβ1Feα1 and His146β2Feβ2Feα2, respectively. Panel B: temporal change of (from top to bottom) the Cα–Cα distance between His146β1 and His146β2, the Cα–Cα distance between His143β1 and His143β2, 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 T0 (2DN2) and oxy R4 (2DN3) states, respectively. Top panel is also given in the main MS.

https://doi.org/10.7554/eLife.35560.004
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 β1, and right panels those for β2. 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.

https://doi.org/10.7554/eLife.35560.005
Figure 1—figure supplement 4
Local structural changes around His146.

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 1μs 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.

https://doi.org/10.7554/eLife.35560.006
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 β1 and β2.

https://doi.org/10.7554/eLife.35560.007
Figure 1—figure supplement 6
Temporal structural changes of the Cα–Cα distance between His146β1 and His146β2 in the 75 Å box for different simulation conditions.

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

https://doi.org/10.7554/eLife.35560.008
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
Figure 3 with 1 supplement
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
Figure 3—figure supplement 1
Number N of H-bonds and their fluctuations δN for pure water and simulations including the protein for the four box sizes.

Left column: Probability distribution of number of H-bonds N relative to the average of the 150 Å box. Right column: probability distribution of the fluctuation δN for each box individually. The left column illustrates the loss in the average number of N relative to the 150 Å box for smaller box sizes. This loss increases with including medium (d=3.0 Å) and weak (d=3.3 Å) 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 δN depends sensitively on the box size.

https://doi.org/10.7554/eLife.35560.012
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
Figure 5 with 4 supplements
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
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 r+R from spherical cavities. R is the (spherical) size of the solute and r 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 25 Å 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 rcoordinate 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 20 Å away from the solute’s surface.

https://doi.org/10.7554/eLife.35560.015
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 g(r) (solid lines) and N(r) (dashed lines) for 90 and 150 Å boxes. The left panel shows the entire range r of distances from the center of the protein and the right panel only that concerning the central cylinder. The value of N(r=8.5) is 150 which is consistent with the results from explicit counting and confirms that the g(r) is meaningful.

https://doi.org/10.7554/eLife.35560.016
Figure 5—figure supplement 3
Distribution P(N) 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.

https://doi.org/10.7554/eLife.35560.017
Figure 5—figure supplement 4
Water at the (α1β2):(α2β1) interface.

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

https://doi.org/10.7554/eLife.35560.018
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

Tables

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

Data availability

All data generated or analysed during this study are included in the manuscript and supporting files.

Additional files

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)