Response to comment on 'Valid molecular dynamics simulations of human hemoglobin require a surprisingly large box size'

Abstract

We recently reported that molecular dynamics simulations for hemoglobin require a surprisingly large box size to stabilize the T(0) state relative to R(0), as observed in experiments (El Hage et al., 2018). Gapsys and de Groot have commented on this work but do not provide convincing evidence that the conclusions of El Hage et al., 2018 are incorrect. Here we respond to these concerns, argue that our original conclusions remain valid, and raise our own concerns about some of the results reported in the comment by Gapsys and de Groot that require clarification.

Keywords: box size; computational biology; hemoglobin; hydrophobic effect; molecular biophysics; molecular dynamics; none; structural biology; systems biology.

Figure 1.. Temporal change of the Cα–Cα distance between His₁₄₆B2 and His1₄₆B2 for the 150 Å box during 1280 nanoseconds of molecular dynamics simulation.

Figure 2.. Reaction time distribution p(τ) for the decomposition HSO₃Cl -> SO₃+HCl from 20 (black), 500 (red) and 5000 (blue) reactive trajectories.

With 20 trajectories the distribution is far from converged and even with 500 trajectories convergence appears to be incomplete. It is only for 5000 independent events that p(τ) approaches a smooth distribution.

Figure 3.. Temporal change in the 150 Å box of (from top to bottom) the Cα–Cα distance between His146B1 and His146B2, the Cα–Cα distance between His143B1 and His143B2, and the Cα RMSD relative to the 2DN2 X-ray structure.

The green and red lines report the time series reported in our article (simulation starting from the T(0) state) and the time series of a simulation starting from the decayed T(0) state (i.e., an R(0) state structure) from R state in the 120 Å box, 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.

Figure 4.. Temporal change of the Cα–Cα distance between His₁₄₆B1 and His₁₄₆B2 for the 75 Å box (A), 90 Å box (B), 120 Å box (C), and 150 Å box (D).

In (A), (C), and (D) the black line reports the time series published and the red line shows the new time series obtained with the new simulation with the Fe–N bond. In (B) the black trace is from El Hage et al. (2018) the red trace for the bonded Fe–His simulations, the green line from using finite cutoffs for evaluating the electrostatics (cutoff 14 Å), and the orange line from using Kovalevsky protonation states for all histidines.

Figure 5.. Time-averaged radial distribution functions g(r) between the C-terminal (COO) of His146 and water H for box sizes 90 Å (top) and 150 Å (middle) data from Figure 6 of El Hage et al. (2018) and from before and after the transition, obtained when using the Kovalevsky protonation in the 90 Å box (bottom).

Figure 6.. Fe-N_His bond lengths in simulations without (left) and with (right) explicit bond.