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

Evaluating sheep hemoglobins with MD simulations as an animal model for sickle cell disease

Sickle cell disease (SCD) affects millions worldwide, yet there are few therapeutic options. To develop effective treatments, preclinical models that recapitulate human physiology and SCD pathophysiology are needed. SCD arises from a single Glu-to-Val substitution at position 6 in the β subunit of hemoglobin (Hb), promoting Hb polymerization and subsequent disease. Sheep share important physiological and developmental characteristics with humans, including the same developmental pattern of fetal to adult Hb switching. Herein, we investigated whether introducing the SCD mutation into the sheep β-globin locus would recapitulate SCD's complex pathophysiology by generating high quality SWISS-MODEL sheep Hb structures and performing MD simulations of normal/sickle human (huHbA/huHbS) and sheep (shHbB/shHbS) Hb, establishing how accurately shHbS mimics huHbS behavior. shHbS, like huHbS, remained stable with low RMSD, while huHbA and shHbB had higher and fluctuating RMSD. shHbB and shHbS also behaved identically to huHbA and huHbS with respect to β2-Glu6 and β1-Asp73 (β1-Asn72 in sheep) solvent interactions. These data demonstrate that introducing the single SCD-causing Glu-to-Val substitution into sheep β-globin causes alterations consistent with the Hb polymerization that drives RBC sickling, supporting the development of a SCD sheep model to pave the way for alternative cures for this debilitating, globally impactful disease.

Implementation of Telescoping Boxes in Adaptive Steered Molecular Dynamics

Long-time dynamical processes, such as those involving protein unfolding and ligand interactions, can be accelerated and realized through steered molecular dynamics (SMD). The challenge has been the extraction of information from such simulations that generalize for complex nonequilibrium processes. The use of Jarzynski's equality opened the possibility of determining the free energy along the steered coordinate, but sampling over the nonequilibrium trajectories is slow to converge. Adaptive steered molecular dynamics (ASMD) and other related techniques have been introduced to overcome this challenge through the use of stages. Here, we take advantage of these stages to address the numerical cost that arises from the required use of very large solvent boxes. We introduce telescoping box schemes within adaptive steered molecular dynamics (ASMD) in which we adjust the solvent box between stages and thereby vary (and optimize) the required number of solvent molecules. We have benchmarked the method on a relatively long α-helical peptide, Ala₃₀, with respect to the potential of mean force and hydrogen bonds. We show that the use of telescoping boxes introduces little numerical error while significantly reducing the computational cost.

Dynamics of camel and human hemoglobin revealed by molecular simulations

Hemoglobin is one of the most widely studied proteins genetically, biochemically, and structurally. It is an oxygen carrying tetrameric protein that imparts the characteristic red color to blood. Each chain of hemoglobin harbors a heme group embedded in a hydrophobic pocket. Several studies have investigated structural variations present in mammalian hemoglobin and their functional implications. However, camel hemoglobin has not been thoroughly explored, especially from a structural perspective. Importantly, very little is known about how the heme group interacts with hemoglobin under varying conditions of osmolarity and temperature. Several experimental studies have indicated that the tense (T) state is more stable than the relaxed (R) state of hemoglobin under normal physiological conditions. Despite the fact that R state is less stable than the T state, no extensive structural dynamics studies have been performed to investigate global quaternary transitions of R state hemoglobin under normal physiological conditions. To evaluate this, several 500 ns all-atom molecular dynamics simulations were performed to get a deeper understanding of how camel hemoglobin behaves under stress, which it is normally exposed to, when compared to human hemoglobin. Notably, camel hemoglobin was more stable under physiological stress when compared to human hemoglobin. Additionally, when compared to camel hemoglobin, cofactor-binding regions of hemoglobin also exhibited more fluctuations in human hemoglobin under the conditions studied. Several differences were observed between the residues of camel and human hemoglobin that interacted with heme. Importantly, distal residues His58 of α hemoglobin and His63 of β hemoglobin formed more sustained interactions, especially at higher temperatures, in camel hemoglobin. These residues are important for oxygen binding to hemoglobin. Thus, this work provides insights into how camel and human hemoglobin differ in their interactions under stress.

The functional role of the hemoglobin-water interface

The interface between hemoglobin (Hb) and its environment, in particular water, is of great physiological relevance. Here, results from in vitro, in vivo, and computational experiments (molecular dynamics simulations) are summarized and put into perspective. One of the main findings from the computations is that the stability of the deoxy, ligand-free T-state (T₀) can be stabilized relative to the deoxy R-state (R₀) only in sufficiently large simulation boxes for the hydrophobic effect to manifest itself. This effect directly influences protein stability and is operative also under physiological conditions. Furthermore, molecular simulations provide a dynamical interpretation of the Perutz model for Hb function. Results from experiments using higher protein concentrations and realistic cellular environments are also discussed. One of the next great challenges for computational studies, which as we show is likely to be taken up in the near future, is to provide a molecular-level understanding of the dynamics of proteins in such crowded environments.

Molecular Dynamics of Hemoglobin Reveals Structural Alterations and Explains the Interactions Driving Sickle Cell Fibrillation

In sickle cell anemia, deoxyhemoglobin deforms RBCs by forming fibrils inside that disintegrate on oxygenation. We studied 100 ns long all-atom molecular dynamics (MD) for sickle and normal hemoglobin fibril models to understand this process, complemented by multiple 1 μs MD for a single tetramer of sickle and normal hemoglobin in deoxy and oxy states. We find that the presence of hydrophobic residues without a bulky side chain at β-6 in hemoglobin is the reason for the stability of the fibrils. Moreover, the free energy landscapes from MD of hemoglobin starting in the tensed (T) state capture the putative transition from T to relaxed (R) state, associated with oxygen binding. The three conformational wells in the landscapes are characterized by the quaternary changes where one αβ dimer rotates with respect to the other. The conformational changes from the oxygenation of sickle hemoglobin hinder the intermolecular contacts necessary for fibril formation.

Ensembles Are Required to Handle Aleatoric and Parametric Uncertainty in Molecular Dynamics Simulation

Classical molecular dynamics is a computer simulation technique that is in widespread use across many areas of science, from physics and chemistry to materials, biology, and medicine. The method continues to attract criticism due its oft-reported lack of reproducibility which is in part due to a failure to submit it to reliable uncertainty quantification (UQ). Here we show that the uncertainty arises from a combination of (i) the input parameters and (ii) the intrinsic stochasticity of the method controlled by the random seeds. To illustrate the situation, we make a systematic UQ analysis of a widely used molecular dynamics code (NAMD), applied to estimate binding free energy of a ligand-bound to a protein. In particular, we replace the usually fixed input parameters with random variables, systematically distributed about their mean values, and study the resulting distribution of the simulation output. We also perform a sensitivity analysis, which reveals that, out of a total of 175 parameters, just six dominate the variance in the code output. Furthermore, we show that binding energy calculations dampen the input uncertainty, in the sense that the variation around the mean output free energy is less than the variation around the mean of the assumed input distributions, if the output is ensemble-averaged over the random seeds. Without such ensemble averaging, the predicted free energy is five times more uncertain. The distribution of the predicted properties is thus strongly dependent upon the random seed. Owing to this substantial uncertainty, robust statistical measures of uncertainty in molecular dynamics simulation require the use of ensembles in all contexts.

Uncertainty quantification in classical molecular dynamics

Molecular dynamics simulation is now a widespread approach for understanding complex systems on the atomistic scale. It finds applications from physics and chemistry to engineering, life and medical science. In the last decade, the approach has begun to advance from being a computer-based means of rationalizing experimental observations to producing apparently credible predictions for a number of real-world applications within industrial sectors such as advanced materials and drug discovery. However, key aspects concerning the reproducibility of the method have not kept pace with the speed of its uptake in the scientific community. Here, we present a discussion of uncertainty quantification for molecular dynamics simulation designed to endow the method with better error estimates that will enable it to be used to report actionable results. The approach adopted is a standard one in the field of uncertainty quantification, namely using ensemble methods, in which a sufficiently large number of replicas are run concurrently, from which reliable statistics can be extracted. Indeed, because molecular dynamics is intrinsically chaotic, the need to use ensemble methods is fundamental and holds regardless of the duration of the simulations performed. We discuss the approach and illustrate it in a range of applications from materials science to ligand-protein binding free energy estimation. This article is part of the theme issue 'Reliability and reproducibility in computational science: implementing verification, validation and uncertainty quantification in silico'.

Biomolecular Modeling and Simulation: A Prospering Multidisciplinary Field

We reassess progress in the field of biomolecular modeling and simulation, following up on our perspective published in 2011. By reviewing metrics for the field's productivity and providing examples of success, we underscore the productive phase of the field, whose short-term expectations were overestimated and long-term effects underestimated. Such successes include prediction of structures and mechanisms; generation of new insights into biomolecular activity; and thriving collaborations between modeling and experimentation, including experiments driven by modeling. We also discuss the impact of field exercises and web games on the field's progress. Overall, we note tremendous success by the biomolecular modeling community in utilization of computer power; improvement in force fields; and development and application of new algorithms, notably machine learning and artificial intelligence. The combined advances are enhancing the accuracy andscope of modeling and simulation, establishing an exemplary discipline where experiment and theory or simulations are full partners.

Rapid, accurate, precise and reproducible ligand-protein binding free energy prediction

A central quantity of interest in molecular biology and medicine is the free energy of binding of a molecule to a target biomacromolecule. Until recently, the accurate prediction of binding affinity had been widely regarded as out of reach of theoretical methods owing to the lack of reproducibility of the available methods, not to mention their complexity, computational cost and time-consuming procedures. The lack of reproducibility stems primarily from the chaotic nature of classical molecular dynamics (MD) and the associated extreme sensitivity of trajectories to their initial conditions. Here, we review computational approaches for both relative and absolute binding free energy calculations, and illustrate their application to a diverse set of ligands bound to a range of proteins with immediate relevance in a number of medical domains. We focus on ensemble-based methods which are essential in order to compute statistically robust results, including two we have recently developed, namely thermodynamic integration with enhanced sampling and enhanced sampling of MD with an approximation of continuum solvent. Together, these form a set of rapid, accurate, precise and reproducible free energy methods. They can be used in real-world problems such as hit-to-lead and lead optimization stages in drug discovery, and in personalized medicine. These applications show that individual binding affinities equipped with uncertainty quantification may be computed in a few hours on a massive scale given access to suitable high-end computing resources and workflow automation. A high level of accuracy can be achieved using these approaches.

On the importance of statistics in molecular simulations for thermodynamics, kinetics and simulation box size

Computational simulations, akin to wetlab experimentation, are subject to statistical fluctuations. Assessing the magnitude of these fluctuations, that is, assigning uncertainties to the computed results, is of critical importance to drawing statistically reliable conclusions. Here, we use a simulation box size as an independent variable, to demonstrate how crucial it is to gather sufficient amounts of data before drawing any conclusions about the potential thermodynamic and kinetic effects. In various systems, ranging from solvation free energies to protein conformational transition rates, we showcase how the proposed simulation box size effect disappears with increased sampling. This indicates that, if at all, the simulation box size only minimally affects both the thermodynamics and kinetics of the type of biomolecular systems presented in this work.

Water Dynamics Around Proteins: T- and R-States of Hemoglobin and Melittin

The water dynamics, as characterized by the local hydrophobicity (LH), is investigated for tetrameric hemoglobin (Hb) and dimeric melittin. For the T₀ to R₀ transition in Hb, it is found that LH provides additional molecular-level insight into the Perutz mechanism, i.e., the breaking and formation of salt bridges at the α₁/β₂ and α₂/β₁ interface is accompanied by changes in LH. For Hb in cubic water boxes with 90 and 120 Å edge length it is observed that following a decrease in LH as a consequence of reduced water density or change of water orientation at the protein/water interface the α/β interfaces are destabilized; this is a hallmark of the Perutz stereochemical model for the T to R transition in Hb. The present work thus provides a dynamical view of the classical structural model relevant to the molecular foundations of Hb function. For dimeric melittin, earlier results by Cheng and Rossky [ Nature 1998, 392, 696-699] are confirmed and interpreted on the basis of LH from simulations in which the protein structure is frozen. For the flexible melittin dimer, the changes in the local hydration can be as much as 30% greater than for the rigid dimer, reflecting the fact that protein and water dynamics are coupled.