Figures and data in Valid molecular dynamics simulations of human hemoglobin require a surprisingly large box size

Figures
Tables
Additional files

6 figures, 1 table and 1 additional file

Figures

Figure 1 with 6 supplements

Download asset Open asset

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 T $_{0}$ (2DN2) and oxy R $_{4}$ (2DN3) states, respectively.

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

Figure 1—figure supplement 1

Download asset Open asset

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

Download asset Open asset

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 $_{β 1}$ Fe $_{β 1}$ Fe $_{α 1}$ and His146 $_{β 2}$ Fe $_{β 2}$ Fe $_{α 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 T $_{0}$ (2DN2) and oxy R $_{4}$ (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

Download asset Open asset

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

Download asset Open asset

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

Download asset Open asset

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

Download asset Open asset

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

Figure 2

Download asset Open asset

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

Download asset Open asset

The average number of H-bonds per water molecule $< N_{H - bonds} >$ /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 $< N_{H - bonds} >$ /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 $< N_{H - bonds} >$ /molecule is insignificant but clearly increases when weak H-bonds are included. We also observe a significant decrease in the fluctuation of $< N_{H - bonds} >$ /molecule between simulations in the smallest and the largest box sizes and a pronounced drop in $< N_{H - bonds} >$ /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

Download asset Open asset

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

Figure 4

Download asset Open asset

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

Download asset Open asset

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

Download asset Open asset

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 $\sim 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 $r -$ 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 $\sim 20$ Å away from the solute’s surface.

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

Figure 5—figure supplement 2

Download asset Open asset

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

Download asset Open asset

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

Download asset Open asset

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

Figure 6

Download asset Open asset

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 (T $_{0}$ ) and 2DN3 (R $_{4}$ ) 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

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

Additional files

Transparent reporting form: https://doi.org/10.7554/eLife.35560.020
Download elife-35560-transrepform-v1.pdf

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)

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

Mendeley

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

Krystel El Hage
Florent Hédin
Prashant K Gupta
Markus Meuwly
Martin Karplus

(2018)

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

eLife 7:e35560.

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

Share this article

Cite this article

Global structural changes depending on box size.

Tetrameric hemoglobin solvated in boxes of different sizes.

Global structural changes depending on box size.

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

Local structural changes around His146.

Structural transition details.

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

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

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

System-size dependence of the water diffusion coefficient.

System-size dependence of the average equilibrium water density.

Enhanced version of Figure 5 from the main MS.

Number of water molecules in the central cylinder of Hb.

Distribution P(N) of water molecules in the central cylinder for the different box sizes.

Water at the (α1β2):(α2β1) interface.

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

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.

Transparent reporting form

Downloads (link to download the article as PDF)

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

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

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

Temporal structural changes of the C $α$ –C $α$ distance between His146 $_{β 1}$ and His146 $_{β 2}$ in the 75 Å box for different simulation conditions.

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

Distribution $P (N)$ of water molecules in the central cylinder for the different box sizes.

Water at the $(α 1 β 2) : (α 2 β 1)$ interface.

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

C $α$ RMSD (in Å) relative to the 2DN2 (T $_{0}$ ) and 2DN3 (R $_{4}$ ) X-ray structure (Park et al., 2006) of the end point structures at 1 $μ$ s.