1. Computational and Systems Biology
  2. Structural Biology and Molecular Biophysics
Download icon

Response to comment on '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, United States
  3. ISIS, Université de Strasbourg, France
  • Cited 1
  • Views 851
  • Annotations
Cite this article as: eLife 2019;8:e45318 doi: 10.7554/eLife.45318

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.

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

Introduction

Our investigation of the box size dependence of hemoglobin in solution was initiated by the disagreement between the experimentally observed stabilization of the unliganded T state relative to an R-like structure (Edelstein, 1971), and molecular dynamics simulations which found that the unliganded T state (T(0)) was not stable and made a transition to an R-like state in times on the order of a hundred nanoseconds (Hub et al., 2010; Yusuff et al., 2012). We began research to determine what could be wrong with the published simulations and found, after looking at many possibilities, that only a surprisingly large (150 Å) simulation box was able to stabilize the T(0) over at least ~1.3 microseconds (El Hage et al., 2018; Figure 1).

Temporal change of the Cα–Cα distance between His146B2 and His146B2 for the 150 Å box during 1280 nanoseconds of molecular dynamics simulation.
https://doi.org/10.7554/eLife.45318.002

Results

Statistics

To estimate rates, a sufficient number of samples is required. An explicit example, although not related to hemoglobin but illustrating the effect of insufficient statistics, is to consider reactive MD simulations for the vibrationally induced photodissociation of HSO3Cl to SO3 and HCl (Figure 2). It is evident that with 20 events the distribution of reaction times is far from converged and does not even allow one to guess the true distribution: see Yosa Reyes et al., 2016 for details on the simulations. This is also why from ab initio MD simulations of vibrationally induced photodissociation, no rates were determined because widely different reaction times for individual trajectories were found when only considering ~100 trajectories (see Figure 2 of Miller and Gerber, 2006). It also remains unclear how the probability of 0.0026 was determined in Gapsys and de Groot (2019) (hereafter referred to as the comment).

Reaction time distribution p(τ) for the decomposition HSO3Cl -> SO3+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.

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

Hydrophobic effect

As stated in the abstract of the El Hage et al. (2018): ‘The results suggest that such a large box is required for the hydrophobic effect, which stabilizes the T(0) tetramer, to be manifested.’ We felt we did not have conclusive evidence for the role of the hydrophobic effect, and have been working with Adam Willard to investigate this point. We note that the comment includes a list of some differences in the simulation setups: these and other differences in the two sets of setups need to be investigated as possible contributors to the difference in the simulation results. A possible problem with the comparisons in Figure 3 of the comment concerns the fact that we only exchanged input files for the 90 Å box; the authors of the comment constructed the larger boxes used in their simulations. In that regard, it is interesting to note that their observed life time for the 90 Å box with its error bars essentially includes the one observed in our simulation. Point on which we do not agree with the constructs used in the comment to reanalyze our results include the following:

# The subvolume analysis which attempts to compare the diffusivity in large boxes but restricted to a smaller subvolume needs to introduce ‘imaginary boundary conditions’ and it is unclear how this was done. Such an approach appears to assume that water outside this artificial volume behaves as bulk water which is not evident. Alternatively, very short analysis times (~2 picoseconds) need to be used to avoid significant displacement of the water molecules (Persson et al., 2017).

# The water-water hydrogen bond network reconfigures when the structure of the protein changes (T vs. R) and from one box size to another box size (see Figure 3 in El Hage et al., 2018). Following a transition in the protein structure these changes occur on the few nanosecond time scale, as a direct response to the protein structure change. This suggests that the protein and water dynamics are coupled and can not be viewed independently.

# As an important test we have solvated the decayed (i.e., R-state) structure from the simulation in the 120 Å box in the 150 Å box and observed that the T-state is partially recovered (Figure 3). Apparently, solvation of the decayed structure in the larger water box provides the driving force required to partially recover the T(0) structure. This is an additional indication that, thermodynamically, T(0) is favored over R(0) in the largest water box without, however, providing information concerning a quantitative measure of the underlying energetics. Full recovery would require much longer simulations.

# Technically, it is desirable to carry out simulations in a setup that allows the solvent – here water – to behave as a pure liquid for quantities such as g(r) or the average number of H-bonds. With hemoglobin as the solute we find this to be the case for a box size of 150 Å.

The normalization procedure of the radial distribution functions used in the comment remains unclear. It should be noted that the radial distribution functions in El Hage et al. (2018) follow the expectations for spherical cavities (Chandler, 2005). Figure 5 in El Hage et al. (2018) aims at highlighting that in a box not sufficiently large, water barely reaches bulk properties as judged from the g(r), whereas it readily does so when using a sufficiently large box (e.g. 150 Å). As no hydrophobic effect is expected for ubiquitin it is not clear what the data in Figure 2—figure supplement 1 of the comment contribute to relating the thermodynamic stability of the T(0) state of hemoglobin with the hydrophobic effect.

Kinetics as a function of box size

This section raises questions as to the statistical significance of the results. As mentioned above and in Figure 2, converged reaction time distributions require hundreds to thousands of transitions. Given the fact that the T(0) state is unstable in the 90 Å and 120 Å boxes it is reasonable that starting from an ensemble of initial structures a distribution of decay times is found. Comparing two such distributions to provide evidence for the presence or absence of a box size dependence requires them to be converged to some degree which probably can not be achieved from 10 to 20 simulations.

It is stated in the comment that with the setup used in our simulations, 20 simulations were done for 90 Å and 150 Å boxes and 10 for a 120 Å box: however, the results shown in Figure 3 of the comment (that the transition is shorter in the larger boxes) is the inverse of what might be expected. This surprising result could be due to poor statistics, to some inconsistency in the use of our setup, or to something we do not understand. Interestingly, the T(0) lifetime in the 90 Å box, for which we had made input files available, is in the range we found (see Figure 1B in El Hage et al., 2018), whereas this is not the case for the 120 Å and 150 Å boxes (for which we did not send input files). We also note a difference in the magnitude of the error bars in Figure 3A of the comment for the different protonation states which requires an explanation. Additional exchanges of information and detailed analyses would serve to clarify these points.

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.

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

When we were still collaborating to investigate the origin of some of the results reported in the comment, we exchanged input files they used in their work and our input files. We note that a previous publication by one of the authors of the comment (Hub et al., 2010) does not provide details on the protonation states of the His, except for His146 for which they used either doubly protonated or protonation at ε. In another publication by one of the authors of the comment (Vesper and de Groot, 2013) it appears that the Kovalevsky protonation states that they sent us were not used. Thus we need more information to verify which protonation states were used in the various simulations. Our comparison showed that the protonation states of the histidine residues were very different (Table 1). We note that it is not sufficient to characterize the histidine residues as ‘protonated’ or ‘deprotonated’. Rather, it is also relevant whether protonation is at the δ or ε positions, for example His143 which should be ε-protonated as it anchors the water network stabilizing His146 (see Figure 6 in El Hage et al., 2018).

Table 1
Histidine protonation states in hemoglobin as used in simulations by the present authors (Zheng et al., 2013; El Hage et al., 2018), and the histidine protonation states in the files that were supplied to us by the authors of the comment (originally from Kovalevsky et al., 2010).
https://doi.org/10.7554/eLife.45318.005
Zheng et al., 2013; El Hage et al., 2018Kovalevsky et al., 2010Zheng et al., 2013; El Hage et al., 2018Kovalevsky et al., 2010
Res.Chain A/CChain AChain CRes.Chain B/DChain BChain D
20HSEHSPHSE2HSEHSEHSE
45HSEHSEHSE63HSEHSPHSE
50HSDHSDHSP77HSEHSEHSE
58HSEHSPHSE92HSD-FeHSD-FeHSD-Fe
72HSDHSPHSP97HSEHSPHSP
87HSD-FeHSD -FeHSD-Fe116HSEHSPHSP
89HSEHSEHSP117HSEHSEHSE
103HSEHSPHSP143HSEHSEHSP
112HSEHSPHSP146HSPHSPHSP
122HSEHSEHSE

Further differences between the simulations concern the shape of the simulation box (cubic in El Hage et al., 2018 vs. dodecahedral), the center of mass of the protein was constrained to avoid periodic boundary effects (El Hage et al., 2018), and the terms for the dihedral parameters for the Fe-proximal histidine differed. Removing the angular center of mass motion is only warned against if the motion of the center of mass itself is not controlled (i.e., if the solute can cross the boundaries of the periodic box). However, this is not the case in our simulations, but appeared to happen in some of those performed by the authors of the comment. Concerning the effect of the box shape on the proteins simulated under periodic boundary conditions, it has been found that the box type can have a statistically significant effect on the outcome of a simulation, and that the magnitude of the effect depends on the protein considered (Wassenaar and Mark, 2006). Although all of the latter may contribute to the difference in the results, we believe that a possible difference in the His protonation states used in the simulations would play a major role and needs to be examined.

Kinetics rather than thermodynamics

A necessary but not sufficient criterion for thermodynamic stability is that the structure considered (here T(0)) does not decay on the time scale of the simulations (e.g., to R(0)). However, transition times from T to R can not be determined from such an approach. For this, either a very large number of independent simulations (on the order of hundreds to thousands; see Figure 2) are required, or free energy simulations need to be carried out. As pointed out in El Hage et al. (2018) from the available experimental data, the T(0) to R(0) transition time should be on the order of seconds. The transitions in simulation boxes of 75 Å, 90 Å, and 120 Å all were under one microsecond and increased with increasing box size (Figure 4). Panels (A), (B), and (D) show that the presence or absence of the Fe-His bond does not lead to substantial changes of the time series as far as they have been followed. Panel (B) also reports the changes for simulations with electrostatic cutoffs (green) and when simulations with the Kovalevsky protonation states were carried out (orange). With finite cutoff the T(0) state is stable for at least 150 nanoseconds, and with the Kovalevsky protonation it decays after 60 nanoseconds, in accord with earlier findings (Hub et al., 2010). The radial distribution functions around His146 in our original study (see Figure 6 in El Hage et al., 2018) and from our simulations using the Kovalevsky protonation in the 90 Å box suggest that with the latter protonation state His146 is typically overhydrated as compared with our (stable) simulation in the 150 Å box (see middle panel in Figure 5). Moreover, with the Kovalevsky protonation in the 90Å box, the height of the first and second peak of g(r) is comparable to that of our own simulation in the 90 Å box after T(0) has decayed (at 500 nanoseconds and later; see Figure 4B).

Temporal change of the Cα–Cα distance between His146B1 and His146B2 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.

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

Comparisons of simulations with and without a Fe–His bond show that the Fe–His bond lengths are comparable in all four subunits (see Figure 6). As to the Fe–Fe distances between the various heme groups, only the 120 Å and 150 Å box simulations gave results in agreement with experiment; the smaller boxes did not.

Fe-NHis bond lengths in simulations without (left) and with (right) explicit bond.
https://doi.org/10.7554/eLife.45318.008

The 150 Å box in a T(0) state was also stable for at least 1.3 microseconds, the length of the simulation (Figure 1).

An obvious improvement of our initial study is to run multiple simulations for each box, in particular if one is interested in the kinetics of the T-to-R transition. We are doing additional simulations to examine the box size dependence of the lifetime. However, we point out that it is very difficult to do a sufficient number of independent simulations to obtain converged statistical results for a system the size of hemoglobin. In ten simulations for the 120 Å box reported in the comment, 4 out of the 10 lasted longer than one microsecond, compared with 3 out of 10 for the 90 Å box. In simulation of kinetics hundreds to thousands of individual, reactive trajectories are required for converged reaction time distributions that can be compared with experiment (Yosa and Meuwly, 2011; Soloviov et al., 2016). The data in Figure 3 of the comment cannot report on the kinetics because the number of events observed is too small (Figure 2).

Kinetics as a function of box size

We refer readers to our author response. We also note that previous work (Wassenaar and Mark, 2006) did find a dependence on box shape.

Thermodynamics as a function of box size

The results for the other two systems are interesting but there is no reason to expect a box size dependence for the structural changes that were studied. Among other points, the surface area change is small and a hydrophobic effect is unlikely. There is no direct relationship to our findings for hemoglobin.

Conclusions

We believe we have demonstrated that the details of our simulations provide evidence that there is a box size dependence in hemoglobin simulations and that it is likely, though not proved, that in a 150 Å box the unliganded T state is stable. The role of the particular choice of His protonation states in the fast decay presented in the comment remains to be evaluated. Finally, we reiterate that the unliganded T state is known to be stable from experiment: therefore, if our results were not correct, there would have to be a yet unknown problem with simulations of unliganded T state hemoglobin.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15

Decision letter

  1. Yibing Shan
    Reviewing Editor; DE Shaw Research, United States
  2. John Kuriyan
    Senior Editor; University of California, Berkeley, United States
  3. Wei Yang
    Reviewer

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.

[Editors’ note: the original decision letter was sent on 20 March 2019, and was followed by further letters on 29 April 2019 and 7 May 2019: the point-by-point response by Gapsys and de Groot was attached to the last of these letters.]

20 March 2019:

Thank you for submitting your response to the comment by Gapsys and de Groot on your eLife article "Valid molecular dynamics simulations of human hemoglobin require a surprisingly large box size". Your response has been reviewed by two peer reviewers, and the evaluation has been overseen by a Reviewing Editor (Yibing Shan), a Senior Editor (John Kuriyan) and the eLife Features Editor (Peter Rodgers).

We plan to progress this matter as follows:

- i) We will ask Gapsys and de Groot to revise their comment to address the points made by the referees (which include a request to delete some material) and to clarify some of the points raised in your response to their comment.

- ii) The comments of the referees on your response are below. No action with respect to these comments is needed at this stage.

- iii) Once we have received the revised version of the comment from Gapsys and de Groot, we will send it to you so that you can revise your response accordingly. At this stage we will also request some editorial revisions and ask for a point-by-point response to the points in the decision letter (as is our standard practice).

29 April 2019:

Further to our email of March 20 (above), we have now received the revised version of the comment by Gapsys and de Groot, so we are writing with the details of the revisions that you need to make to your response to this comment. As indicated in the email of March 20, you will need to make: i) editorial changes; ii) changes in response to the reviewer comments; iii) changes to reflect the changes in the comment by Gapsys and de Groot. Please see below for more details.

Please be aware that if your manuscript is accepted for publication, we will publish the decision letter and your point-by-point response to it (excluding minor comments).

# Editorial changes

a) Please change the title to the following (which is the standard format for such articles in eLife):

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

b) Please revise the abstract so that it better summarizes your manuscript. I would suggest something along the following lines:

"In a previous paper we reported that [QUERY: Please summarize your original paper and the dependence on box size here] (El Hage et al., 2018). Gapsys and de Groot have expressed a number of concerns about this work, and have argued that the dependence on box size that we reported is not reproducible (Gapsys and de Groot, 2019). Here we respond to these concerns and argue that our original conclusions remain valid."

Please feel free to revise the above, but please avoid stating that Gapsys and de Groot are incorrect. Please also answer the query.

c) Please revise the conclusion so that, again, it better summarizes your manuscript. I would suggest something along the following lines:

"We continue to believe that our work provides evidence for a dependence on box size in molecular dynamics simulations of human hemoglobin, and that is likely, though not proven, that the unliganded T state of human hemoglobin is stable in a 150 Å box. We also believe that the role of the His protonation states in the fast decay of the T(0) state [OK?] reported by Gapsys and de Groot remains to be evaluated. Finally, we reiterate that the unliganded T state is known to be stable from experiment: therefore, if our results are not correct, there exists an as yet unknown problem with the simulations of unliganded T state hemoglobin."

Please feel free to revise the above.

# Changes in response to the reviewer comments

The reviewer comments (which were included in the email of 20 March) are copied below. Of the six points made by reviewer #1, points 1 and 5 are covered below (under the heading "Changes to reflect the changes in the comment by Gapsys and de Groot"), and points 2, 3 and 4 do not require a response. The four paragraphs from reviewer #2 do not require a response.

Reviewer #1:

This paper is to reply to the comment by Gapsys and de Groot on the published work "'Valid molecular dynamics simulations of human hemoglobin require a surprisingly large box size' by El Hage et al. The reply was focused on that (a) the examination work does not invalidate the original conclusion that there is a significant box size effect and (b) several statements in the comment need to be clarified. With careful review of both the comment by Gapsys and de Groot and this reply, by and large, I agree with the general discussion of El Hage et al., particularly on the fact that the observed box size dependency is interesting and needs to be further explored; and it is likely an explanation to the mystery of unexpected fast T-R transition in MD simulations. Given that, I have the following opinions and suggestions:

1) I have suggested some major changes to Gapsys and de Groot regarding their focus of the discussion. Hopefully they will revise the comment accordingly and the authors of this reply can change their reply accordingly.

2) I appreciate the fact that El Hage et al. further clarified the fact that the hydrophobic effect was proposed as an hypothesis rather than a firm conclusion and they are further exploring this effect.

3) I appreciate the fact that El Hage et al. further clarified that the focus of the original paper is to "demonstrate that the details of the simulations provide evidence that there is a box size dependence in hemoglobin simulations and that is likely, though not proved, that in a 150 Å box the unliganded T state is stable."

4) In my opinion, the importance of the original work even with the above two clarifications still firmly holds and the paper will be an important literature in the field of MD simulation.

5) I agree on the proposal by El Hage et al. that to re-exam their original observation, the simulation details need to be carefully checked. I would like to encourage the authors of both sides to work together to further explore these important observations in the near future.

6) A minor suggestion is that El Hage et al. may consider putting some appreciation comments on the fact that their work has attracted attentions of the authors of the comment; and Gapsys and de Groot have put quite a bit efforts in further looking to this issue and provided some suggestion on some aspect of the analysis. It is a healthy scientific process.

Reviewer #2:

While it is an interesting and potentially important possibility that simulation box size may affect simulation outcome, I think it is clear that the statistics in the paper by El Hage et al. are insufficient to establish this. In the comment by Gapsys and de Groot, more repeats of the simulation of hemoglobin were reported, but the number of the repeats is not large enough to exclude this possibility, although the simulations appear to be inconsistent with this possibility. The suggested effect of simulation box size is necessarily a subtle one, if it does exist, and the value of the comment, in my assessment, is in pointing out that this is an unresolved issue that require further investigations. This, importantly, is acknowledged by El Hage et al in their reply to the comment.

Gapsys and de Groot correctly pointed out that even the transition time is established to be smaller in simulations of a small box than in ones of a big box, one cannot reliably infer that the thermostability of the initial T state is reduces by the small box size. However, this is common knowledge of even beginners of this field, of which surely El Hage et al. are aware. This distinction between kinetics and thermodynamics is peripheral at this point when it is not known for certain whether the transition time is indeed affected by box size. I think this section should be removed from the comment.

Gapsys and de Groot proposed an alternative approach to analyze the potential hydrophobic effect associated with box size and showed that such an approach would lead to a different conclusion. This is valid. Neutral statements to this effect would suffice.

The suggested effect of simulation box size is necessarily a subtle one. Thus the fact that Gapsys and de Groot did not see it in simulations of a number of other systems is not sufficient to exclude this possibility. It seems unnecessary to report these simulations in the comment. Or at least, it should be acknowledged that these additional simulations are not conclusive either way.

# Changes to reflect changes in the Comment by Gapsys and de Groot

Please address the following points by revising your manuscript and explaining in your point-by-point response how your manuscript has been revised. It is important that you make these changes so that readers of the two articles can get a clear idea of what has been done by both groups.

i) The order of subsections within the Results section of the comment by Gapsys and de Groot has been changed, so please change the order of the subsections within the Results section of your manuscript.

ii) Gapsys and de Groot have expanded the "Statistics" subsection by revising the sentences that start "For example, we find.. . " and "This is further underscored... " They have also added eight new sentences to the end of this section, and have added a new figure (Figure 1-figure supplement 1).

Please consider revising your manuscript to reflect these additions.

iii) Gapsys and de Groot have expanded the "Hydrophobic effect" subsection by adding six new sentences to the end of this sectio. They have also added a new figure (Figure 2-figure supplement 1).

Please consider revising your manuscript to reflect these additions.

iv) Gapsys and de Groot have added a new sentence to the end of the "Kinetics rather than thermodynamics" subsection.

Please consider revising your manuscript to reflect this addition.

v) We asked Gapsys and de Groot to revise the previous version of this manuscript to better explain how their setups compared to/differed from your setups. Our requests to them are in regular type below; their replies to these requests are in italic type. Please revise your manuscript as necessary in the light of these changes:

"We note that the setups by El Hage et al and by Gapsys and de Groot are not identical, and the differences in the setups potentially involve several parameters (e.g., how the initial water box was constructed, how long the water molecules were equilibrated before the production runs, whether the protein system was restrained during the water equilibriation, and what types of baro-stat and thermo-stat were used et al.) To clarify these issues, please be explicit about the differences between your setups and the setups used by El Hage et al.

REPLY: We are aware of this issue and therefore had contacted El Hage et al for their setup. The data in our comment that refer to 'El Hage setup' is based on their 90A setup and therefore is identical to that of El Hage et al in all the mentioned aspects of how the initial water box was constructed, how long the water molecules were equilibrated before the production runs, how the protein system was restrained during the water equilibriation, what force field was applied, and what types of barostat and thermostat et al were used. Therefore, the reviewers' assertion above "We note that the setups by El Hage et al and by Gapsys and de Groot are not identical" appears incorrect. Based on this setup of the 90A box, we ran 20 simulations of 1 microsecond each. 9 of these did not complete the transition, and thus remain in the T state, which in the El Hage et al. 2018 study had been claimed as an exclusive feature of the 150A box. The fact that we observe the lack of a transition also in equally long simulations in the 90A box thus renders the conclusion of El Hage et al. invalid that the lack of a transition to R (or indeed a valid description of deoxyhemoglobin to remain 'stable' in T) is due to the larger box size in the 150A box.

In addition, we have also carried out simulations in larger boxes of cubic sizes of 120A and 150A. The 120A and 150A boxes were constructed from the 90A box by adding a water layer (with solvated ions) of 30 and 60A, respectively. These therefore also are identical in terms of protein configuration as well as protein solvation shell to the El Hage et al setup. We have now added this description to the manuscript. That the added solvation has created a setup indistinguishable from the El Hage et al setup for the 120A and 150A boxes is illustrated by the fact that the single trajectory findings of El Hage et al. fall well within the error bars that we provide not only for the 90A box but also for the 120A and 150A boxes. This can be even more clearly seen in our Fig. 1B, that shows that the El Hage et al. simulations fall well within the distributions we find for all the investigated box sizes.

That other parameters in the setup indeed have an influence on the transition statistics is reflected in Fig. 3A where we show that three different hemoglobin setups in terms of e.g. chosen protonation states or protein-heme interactions show three different types of transition kinetics. In none of these three, we see a systematic effect of the box size, however.

Also, please revise your manuscript to address the following comments in the response from El Hage et al. a) The comments in the sentence that starts "A possible problem with the comparisons in their figure 3…"

REPLY: The setup issue is already addressed above. It was the choice of El Hage et al to only send us their 90A setup (we had before that shared all our setups with them). The only thing we did to construct the 120A and 150A boxes was to add a layer of solvent with ions around the pre-solvated 90A box. In the same paragraph, El Hage et al write: "In that regard, it is interesting to note that their observed life time for the 90A box with its error bars essentially includes the one observed in our simulation." The same is in fact true for the 120A and 150A boxes: the single trajectory findings of El Hage et al. fall well within the error bars that we provide. This can be even more clearly seen in our Fig. 1B, that shows that the El Hage et al. simulations fall well within the distributions we find. There is thus no reason to assume 'A possible problem' with the 120A and 150A boxes.

b) The comments in the sentence that starts "We also note a difference in the magnitude of the error bars.. . "

REPLY: the errors are different because the transition times are different. Smaller absolute transition times result in smaller associated absolute errors. We now included this notion in the revised manuscript.

c) The comments in the paragraph that starts "When we were still collaborating…"

REPLY: This section is about the effect of protonation states. For the systems labeled 'El Hage' in our comment, the protonation states are identical to the ones of El Hage et al (2018) as the simulations started from the input provided by El Hage et al. In addition, as pointed out in our comment, The El Hage et al protonation states are nearly identical to the ones in our "Hub" setup (except for one small difference of a surface histidine of which we argue that it is unlikely to substantially affect the transition kinetics) and therefore it seems odd to state (like El Hage et al do): "protonation states of the histidine residues were very different". They are different to the Kovalevsky et al. setup, the effect of which on hemoglobin we had already reported in 2013 (Vesper and de Groot, 2013). This section in the reply of El Hage et al thus seems rather confusing and does not appear to provide novel insight. In addition, the protonation states listed in Table 1 of the El Hage et al reply differ from the setup that was sent to us.

In total, we tested 3 different hemoglobin setups with different protonation states and observed no box size dependence in any of them."

7 May 2019:

Thank you for your email. I attach the point-by-point response from Gapsys and de Groot in response to the referee reports on the previous version of their manuscript. Please treat this document in confidence (and please do not quote from it in your revised manuscript).

As requested in my email of April 29, can you please revise your manuscript to i) address the editorial points outlined in that letter; ii) address the comments from the referees included in that letter; iii) to reflect the changes in the manuscript of Gapsys and De Groot. Please note that it should be possible to address all these points by making textual revisions to your manuscript: that is, we are not asking you to perform further simulation.

You will be able to use the following link to upload your revised manuscript and your point-by-point response to the email of April 29:

I look forward to receiving your revised manuscript and point-by-point response. It would be helpful if you could submit these within the next two weeks (i.e., by Tuesday 21 May). If you think you will need more time to revise your manuscript, please let me know.

Finally, please note that eLife publishes the decision letter and the point-by-point response it (excluding minor comments) for all papers.

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

Author response

In response to the comments in the decision letter we have made the following changes to our Response to the Comment from Gapsys and de Groot.

# Editorial changes

a) Please change the title to the following (which is the standard format for such articles in eLife): Response to comment on ‘Valid molecular dynamics simulations of human aemoglobin require a surprisingly large box size’

REPLY: This was done.

b) Please revise the abstract so that it better summarizes your manuscript. I would suggest something along the following lines: “In a previous paper we reported that [QUERY: Please summarize your original paper and the dependence on box size here] (El Hage et al., 2018). Gapsys and de Groot have expressed a number of concerns about this work, and have argued that the dependence on box size that we reported is not reproducible (Gapsys and de Groot, 2019). Here we respond to these concerns and argue that our original conclusions remain valid." Please feel free to revise the above, but please avoid stating that Gapsys and de Groot are incorrect. Please also answer the query.

REPLY: The abstract was partly rewritten. We did not state that Gapsys and de Groot are incorrect, hence, no correction was required.

c) Please revise the conclusion so that, again, it better summarizes your manuscript. I would suggest something along the following lines: “We continue to believe that our work provides evidence for a dependence on box size in molecular dynamics simulations of human hemoglobin, and that is likely, though not proven, that the unliganded T state of human hemoglobin is stable in a 150 box. We also believe that the role of the His protonation states in the fast decay of the T(0) state [OK?] reported by Gapsys and de Groot remains to be evaluated. Finally, we reiterate that the unliganded T state is known to be stable from experiment: therefore, if our results are not correct, there exists an as yet unknown problem with the simulations of unliganded T state hemoglobin.” Please feel free to revise the above.

REPLY: The Conclusion was partly rewritten.

# Changes in response to the reviewer comments

Reviewer 1:

1) I have suggested some major changes to Gapsys and de Groot regarding their focus of the discussion. Hopefully they will revise the comment accordingly and the authors of this reply can change their reply accordingly.

REPLY: Apparently the reviewer had concerns about the Comment by Gapsys and de Groot, which were not addressed since no changes were made by Gapsys and de Groot in their Comment, other than to report additional simulations.

REPLY: Please see below for our comments on the point-by-point response from Gapsys and de Groot.

5) I agree on the proposal by El Hage et al. that to re-exam their original observation, the simulation details need to be carefully checked. I would like to encourage the authors of both sides to work together to further explore these important observations in the near future.

REPLY: Please see for our comments on the point-by-point response from Gapsys and de Groot (below).

Reviewer 2: no response required

# Changes to reflect changes in the Comment by Gapsys and de Groot

i) The order of subsections within the Results section of the comment by Gapsys and de Groot has been changed, so please change the order of the subsections within the Results section of your manuscript.

REPLY: This was done.

ii) Gapsys and de Groot have expanded the “Statistics” subsection by revising the sentences that start “For example, we find.. . ” and “This is further underscored... “. They have also added eight new sentences to the end of this section, and have added a new figure (Figure 1-figure supplement 1). Please consider revising your manuscript to reflect these additions.

REPLY: We expanded our own comments on this point and provide quantitative information that from 10 or 20 reactive events the kinetics of a process cannot be reliably described.

iii) Gapsys and de Groot have expanded the “Hydrophobic effect” subsection by adding six new sentences to the end of this section. They have also added a new figure (Figure 2-figure supplement 1). Please consider revising your manuscript to reflect these additions.

REPLY: We included our reservations towards the normalization procedure which is still not satisfactorily explained. We also comment directly on the new Figure 2-figure supplement 1.

iv) Gapsys and de Groot have added a new sentence to the end of the “Kinetics rather than thermodynamics” subsection. Please consider revising your manuscript to reflect this addition.

REPLY: As explained in the 'Statistics' subsection, the number of events (10 to 20) is too small to provide reliable information about the kinetics, see also our new Figure 2. We added a comment on that to the 'Kinetics rather than thermodynamics' subsection.

v) We asked Gapsys and de Groot to revise the previous version of this manuscript to better explain how their setups compared to/differed from your setups. Our requests to them are in regular type below; their replies to these requests are in italic type. Please revise your manuscript as necessary in the light of these changes.

REPLY: No change was made in our Response to the comment. See also for our comments on the point-by-point response from Gapsys and de Groot (below)

We thus have addressed all points raised by the reviewers and we hope that the work can now be accepted for publication.

[Editors’ note: in the rest of this document, Profs Muewly and Karplus respond to the author response from Gapsys and de Groot.]

Point 1.

In point 1, the reviewer focuses on the possibility that the box size effect on the conformational transition of large proteins cannot be excluded, although it is far from established by the work of El Hage et al, 2018. The Comment from Gapsys and de Groot reports more repeats of simulations (approximately 10), than El Hage et al. (just one), but it also falls short of establishing the statistical significance of an opposite observation for which hundreds or more simulations would be required. We have provided concrete data for a much smaller system in Figure 2 in our Response to the Comment (El Hage et al., 2019).

"The reviewers believe it is important for both parties to not overstate the case, and to clearly acknowledge that the box size effect remains a hypothesis, neither established nor convincingly precluded. Unless the authors think this is not the case, the comment should be revised to reflect this basic reality that neither parties in this debate knows the ultimate answer here."

In response to this comment of the reviewers, Gapsys and de Groot introduce a comment about the very low probability of longer lifetimes but present no justification for this statement. Also, they quote a paper, which claims that only 10 repeats are required to obtain a statistically significant result. This contrasts sharply with several published papers (Yosa and Meuwly, 2011; Yosa Reyes et al., 2016; see also Figure 2 in El Hage et al., 2019) where thousands of events were found to be required to describe the kinetics. For a comparison, see Figure 2 in Miller and Gerber, 2006, where 100 simulations were not sufficient to converge the reaction probability.

Point 2.

The reviewer states "The discussion on the distinction between thermal and kinetic stability in the first paragraph of the results section is not needed and should be removed. (The corresponding section in the reply by El Hage et al. should be also be removed.)"

Here we agree with Gapsys and de Groot that the question of whether a kinetic or thermodynamic effect is involved is important and should be discussed. However, this should only be addressed once the simulation conditions have been established that lead to a thermodynamically stable T(0) state, consistent with experiment. Moreover, we disagree with the rationale of Gapsys and de Groot for doing so. The purpose of El Hage et al. (2018) was never to quantify the actual kinetics (i.e. transition time). In fact, it is stated in El Hage et al., 2018 that "..the R(0) to T(0) transition time is about 20 $\mu$s, much longer than the simulation time."

El Hage et al., 2018 reported on a possible simulation box size effect on the T(0) stability and was motivated, as mentioned at the beginning of this reply, by several published simulations of the T(0) state of hemoglobin, of which Hub et al., 2010) was one. As stated in the Introduction of El Hage et al., 2018, it is known from experiment that T(0) is stable so that the observation in the simulations of Hub et al., 2010 that the T(0) state goes over to an R(0)-like state indicates that there is something wrong with the thermodynamics. Further, that the T(0) to R(0) transition takes place on the hundred nanosecond time scale, while it is known from experiment (see El Hage et al., 2018) that such rare transitions require seconds, indicates that something is also wrong with the kinetics in the simulations of Gapsys and de Groot.

Point 3.

The reviewer begins with the statement "It is known though rarely formally acknowledged, that when the simulation length is far shorter than the actual time scale of a molecular event, the simulation results are sensitive to the exact simulation set up. We note that the setup of El Hage and by Gapsys and de Groot are not identical and the differences in the setup potentially involve several parameters."

The reviewer then lists several differences between the simulation setups.

In their response, Gapsys and de Groot point out that for the 90Å box simulated with Gromacs we had sent them our input files (as we believed we were still collaborating), and so it is likely that the simulation conditions were identical. However, for the 120Å and 150Å boxes they did not have our input files and so they generated the larger boxes in a rather arbitrary manner by "adding a water layer (with solvated ions of 30 Angstroms and 60 Angstroms, respectively." There is no reason to presume that this way of creating the larger boxes yielded a setup that is identical to the ones we used. Given the coupling between the protein conformation and water (pointed out on page 4 of the Reply to the Comment and elaborated under point 4), it is very likely that the setups were, in fact, significantly different.

We also note here that Gapsys and de Groot list a set of differences between our simulations and theirs that could be important since the effect is subtle. Reviewer 1 suggests further communication to resolve the difference. Gapsys and de Groot declined to do so.

As to specific points of the reviewer where Gapsys and de Groot are asked to revise their manuscript

- point (a) See comments above about the 120Å and 150Å boxes.

- point (b) the response of Gapsys and de Groot is appropriate.

- point (c) See the comments about the protonation states used in the simulations in the section 'Kinetics as a function of box size' in the main article.

Point 4.

This comment refers to the hydrophobic effect and the reviewers "suggest that in the revision the authors make clear that the hydrophobic effect hypothesis is not precluded although the original analysis by El Hage can be improved." Gapsys and de Groot argue that the various effects observed by El Hage, which suggested that there is a box size dependence, can be explained as a simple dilution effect. However, as stated in our Reply, "we do not agree with the constructs used in the Comment to reanalyze our results. These include: The subvolume analysis which attempts to compare the diffusivity in large boxes but restricted to a smaller subvolume needs to introduce 'imaginary boundary conditions' and it is unclear how this was done. Such an approach also appears to assume that water outside this artificial volume behaves as bulk water. Alternatively, very short analysis times (~2 ps) need to be used to avoid significant displacement of the water molecules (Persson et al., 2017)."

The water-water hydrogen bond network reconfigures when the structure of the protein changes (T vs. R) and from one box size to another box size (see Figure 3 in El Hage et al., 2018). Following a transition in the protein structure these changes occur on the few nanosecond time scale, as a direct response to the protein structure change. This suggests that the protein and water dynamics are coupled and can not be viewed independently."

This second point is an essential element that was neglected in the constructs of Gapsys and de Groot and makes clear that the observed behavior cannot be ascribed to a simple dilution effect. Gapsys and de Groot responded: "...box size effect on the hydrophobic effect (as quantified by hydrogen bonding, solvent RDF..) disappears when normalized.." This basically says - or at least implies - that the analysis in Chandler, 2005 is incorrect as this work links the hydrophobic effect to the behaviour of the RDF and our analysis for the different box sizes reflect exactly that behaviour of the RDF. So something with the "normalization" in Gapsys and de Groot must be wrong.

Point 5.

The reviewer points out that "The authors' simulations of small proteins without observing a box size effect does not serve as proof to preclude the box size effect, especially for large systems. The authors should make this important caveat clear in their discussion." In fact, reviewer 2 in his reply to us states "… it seems unnecessary to report these simulations in the Comment. Or at least … are not conclusive either way." Gapsys and de Groot state that they limit their findings to "for the investigated systems". Thus, they do mention the limitation, but the emphasis of what they wrote and what the reviewer requests are different.

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

Article and author information

Author details

  1. Krystel El Hage

    Department of Chemistry, University of Basel, Basel, Switzerland
    Contribution
    Investigation, 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
    Writing—review and editing
    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
    Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4792-7538
  4. Markus Meuwly

    Department of Chemistry, University of Basel, Basel, Switzerland
    Contribution
    Conceptualization, Resources, Formal analysis, Supervision, Validation, Methodology, Writing—original draft, 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, United States
    2. Laboratoire de Chimie Biophysique, ISIS, Université de Strasbourg, Strasbourg, France
    Contribution
    Writing—original draft, Writing—review and editing
    For correspondence
    marci@tammy.harvard.edu
    Competing interests
    No competing interests declared

Funding

Swiss National Science Foundation (200021-117810)

  • Markus Meuwly

Swiss National Science Foundation (NCCR MUST)

  • Markus Meuwly

CHARMM Developent Project

  • Martin Karplus

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

Senior Editor

  1. John Kuriyan, University of California, Berkeley, United States

Reviewing Editor

  1. Yibing Shan, DE Shaw Research, United States

Reviewer

  1. Wei Yang

Publication history

  1. Received: January 21, 2019
  2. Accepted: June 6, 2019
  3. Version of Record published: June 20, 2019 (version 1)

Copyright

© 2019, 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

  • 851
    Page views
  • 45
    Downloads
  • 1
    Citations

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

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. Computational and Systems Biology
    2. Stem Cells and Regenerative Medicine
    Alexander J Tarashansky et al.
    Tools and Resources
    1. Computational and Systems Biology
    2. Microbiology and Infectious Disease
    Michael R McLaren et al.
    Research Article