Molecular dynamics study of the photodissociation of carbon monoxide from myoglobin: Ligand dynamics in the first 10 ps

CHARMM at 45: Enhancements in Accessibility, Functionality, and Speed

Since its inception nearly a half century ago, CHARMM has been playing a central role in computational biochemistry and biophysics. Commensurate with the developments in experimental research and advances in computer hardware, the range of methods and applicability of CHARMM have also grown. This review summarizes major developments that occurred after 2009 when the last review of CHARMM was published. They include the following: new faster simulation engines, accessible user interfaces for convenient workflows, and a vast array of simulation and analysis methods that encompass quantum mechanical, atomistic, and coarse-grained levels, as well as extensive coverage of force fields. In addition to providing the current snapshot of the CHARMM development, this review may serve as a starting point for exploring relevant theories and computational methods for tackling contemporary and emerging problems in biomolecular systems. CHARMM is freely available for academic and nonprofit research at https://academiccharmm.org/program.

Binding of Carbon Monoxide to Hemoglobin in an Oxygen Environment: Force Field Development for Molecular Dynamics

Carbon monoxide (CO) is a byproduct of the incomplete combustion of carbon-based fuels, such as wood, coal, gasoline, or natural gas. As incomplete combustion in a fire accident or in an engine, massively produced CO leads to a serious life threat because CO competes with oxygen (O2) binding to hemoglobin and makes people suffer from hypoxia. Although there is hyperbaric O2 therapy for patients with CO poisoning, the nanoscale mechanism of CO dissociation in the O2-rich environment is not completely understood. In this study, we construct the classical force field parameters compatible with the CHARMM for simulating the coordination interactions between hemoglobin, CO, and O2, and use the force field to reveal the impact of O2 on the binding strength between hemoglobin and CO. Density functional theory and Car-Parrinello molecular dynamics simulations are used to obtain the bond energy and equilibrium geometry, and we used machine learning enabled via a feedforward neural network model to obtain the classical force field parameters. We used steered molecular dynamics simulations with a force field to characterize the mechanical strength of the hemoglobin-CO bond before rupture under different simulated O2-rich environments. The results show that as O2 approaches the Fe2+ of heme at a distance smaller than ∼2.8 Å, the coordination bond between CO and Fe2+ is reduced to 50% bond strength in terms of the peak force observed in the rupture process. This weakening effect is also shown by the free energy landscape measured by our metadynamics simulation. Our work suggests that the O2-rich environment around the hemoglobin-CO bond effectively weakens the bonding, so that designing of O2 delivery vector to the site is helpful for alleviating CO binding, which may shed light on de novo drug design for CO poisoning.

Characterization of the Bottlenecks and Pathways for Inhibitor Dissociation from [NiFe] Hydrogenase

[NiFe] hydrogenases can act as efficient catalysts for hydrogen oxidation and biofuel production. However, some [NiFe] hydrogenases are inhibited by gas molecules present in the environment, such as O2 and CO. One strategy to engineer [NiFe] hydrogenases and achieve O2- and CO-tolerant enzymes is by introducing point mutations to block the access of inhibitors to the catalytic site. In this work, we characterized the unbinding pathways of CO in the complex with the wild-type and 10 different mutants of [NiFe] hydrogenase from Desulfovibrio fructosovorans using τ-random accelerated molecular dynamics (τRAMD) to enhance the sampling of unbinding events. The ranking provided by the relative residence times computed with τRAMD is in agreement with experiments. Extensive data analysis of the simulations revealed that from the two bottlenecks proposed in previous studies for the transit of gas molecules (residues 74 and 122 and residues 74 and 476), only one of them (residues 74 and 122) effectively modulates diffusion and residence times for CO. We also computed pathway probabilities for the unbinding of CO, O2, and H2 from the wild-type [NiFe] hydrogenase, and we observed that while the most probable pathways are the same, the secondary pathways are different. We propose that introducing mutations to block the most probable paths, in combination with mutations to open the main secondary path used by H2, can be a feasible strategy to achieve CO and O2 resistance in the [NiFe] hydrogenase from Desulfovibrio fructosovorans.

Liquid-Phase Effects on Adsorption Processes in Heterogeneous Catalysis

Aqueous solvation free energies of adsorption have recently been measured for phenol adsorption on Pt(111). Endergonic solvent effects of ∼1 eV suggest solvents dramatically influence a metal catalyst's activity with significant implications for the catalyst design. However, measurements are indirect and involve adsorption isotherm models, which potentially reduces the reliability of the extracted energy values. Computational, implicit solvation models predict exergonic solvation effects for phenol adsorption, failing to agree with measurements even qualitatively. In this study, an explicit, hybrid quantum mechanical/molecular mechanical approach for computing solvation free energies of adsorption is developed, solvation free energies of phenol adsorption are computed, and experimental data for solvation free energies of phenol adsorption are reanalyzed using multiple adsorption isotherm models. Explicit solvation calculations predict an endergonic solvation free energy for phenol adsorption that agrees well with measurements to within the experimental and force field uncertainties. Computed adsorption free energies of solvation of carbon monoxide, ethylene glycol, benzene, and phenol over the (111) facet of Pt and Cu suggest that liquid water destabilizes all adsorbed species, with the largest impact on the largest adsorbates.

Stationary and Time-Dependent Carbon Monoxide Stretching Mode Features in Carboxy Myoglobin: A Theoretical-Computational Reappraisal

The stationary and time-dependent infrared spectrum (IR) of the CO stretching mode (ν_CO) in carboxymyoglobin (MbCO), a longstanding problem of biophysical chemistry, has been modeled through a theoretical-computational method specifically designed for simulating quantum observables in complex atomic-molecular systems and based on a combined application of long time scale molecular dynamics simulations and quantum-chemical calculations. This study is basically focused on two aspects: (i) the origin of the stationary IR substates (termed as A₀, A₁, and A₃) and (ii) the modeling and the interpretation of the ν_CO energy relaxation. The results, strengthened by a more than satisfactory agreement with the experimental data, concisely indicate that (i) the conformational His64-FeCO relevant substates, i.e., characterized by the formation-disruption of the H-bond between the above moieties, are the main responsible of the presence of two distinct and well separated (A₀ and A₁/A₃) spectroscopic regions; (ii) the characteristic bimodal shape of the A₁/A₃ spectral region, according to our model, is the result of the fluctuation of the electric field pattern as provided by the protein-solvent framework perturbing the bound His64-CO-Heme complex; and (iii) the electric field pattern, in conjunction with the relatively high density of MbCO vibrational states, is also the main determinant of the ν_CO energy relaxation, characterizing its kinetic efficiency.

Femtosecond X-ray spectroscopy of haem proteins

We discuss our recently reported femtosecond (fs) X-ray emission spectroscopy results on the ligand dissociation and recombination in nitrosylmyoglobin (MbNO) in the context of previous studies on ferrous haem proteins. We also present a preliminary account of femtosecond X-ray absorption studies on MbNO, pointing to the presence of more than one species formed upon photolysis.

An electrostatic energy-based charge model for molecular dynamics simulation

The interactions of the polar chemical bonds such as C=O and N-H with an external electric field were investigated, and a linear relationship between the QM/MM interaction energies and the electric field along the chemical bond is established in the range of weak to intermediate electrical fields. The linear relationship indicates that the electrostatic interactions of a polar group with its surroundings can be described by a simple model of a dipole with constant moment under the action of an electric field. This relationship is employed to develop a general approach to generating an electrostatic energy-based charge (EEC) model for molecules containing single or multiple polar chemical bonds. Benchmark test studies of this model were carried out for (CH₃)₂-CO and N-methyl acetamide in explicit water, and the result shows that the EEC model gives more accurate electrostatic energies than those given by the widely used charge model based on fitting to the electrostatic potential (ESP) in direct comparison to the energies computed by the QM/MM method. The MD simulations of the electric field at the active site of ketosteroid isomerase based on EEC demonstrated that EEC gave a better representation of the electrostatic interaction in the hydrogen-bonding environment than the Amber14SB force field by comparison with experiment. The current study suggests that EEC should be better suited for molecular dynamics study of molecular systems with polar chemical bonds such as biomolecules than the widely used ESP or RESP (restrained ESP) charge models.

Adsorption of CO and N2 molecules at the surface of solid water. A grand canonical Monte Carlo study

The adsorption of carbon monoxide and nitrogen molecules at the surface of four forms of solid water is investigated by means of grand canonical Monte Carlo simulations. The trapping ability of crystalline Ih and low-density amorphous ices, along with clathrate hydrates of structures I and II, is compared at temperatures relevant for astrophysics. It is shown that when considering a gas phase that contains mixtures of carbon monoxide and nitrogen, the trapping of carbon monoxide is favored with respect to nitrogen at the surface of all solids, irrespective of the temperature. The results of the calculations also indicate that some amounts of molecules can be incorporated in the bulk of the water structures, and the molecular selectivity of the incorporation process is investigated. Again, it is shown that incorporation of carbon monoxide is favored with respect to nitrogen in most of the situations considered here. In addition, the conclusions of the present simulations emphasize the importance of the strength of the interactions between the guest molecules and the water network. They indicate that the accuracy of the corresponding interaction potentials is a key point, especially for simulating clathrate selectivity. This highlights the necessity of having interaction potential models that are transferable to different water environments.

Electrostatic Potential Optimized Molecular Models for Molecular Simulations: CO, CO2, COS, H2S, N2, N2O, and SO2

Molecular simulations have been widely employed in the discovery of nanoporous materials, such as metal-organic frameworks (MOFs) and zeolite, for energy- and environment-related applications. To achieve simulation predictions with better accuracy, we herein present a collection of molecular models, including carbon monoxide (CO), carbon dioxide (CO₂), carbonyl sulfide (COS), hydrogen sulfide (H₂S), nitrogen (N₂), nitrous oxide (N₂O), and sulfur dioxide (SO₂). These models, denoted as electrostatic potential optimized molecular models (ESP-MMs), are systematically developed to not only reproduce experimental vapor-liquid equilibrium but also have accurate electrostatic potential representation surrounding the molecules. Our results show that, with accurate electrostatic potential representations, ESP-MMs can offer improved predictions in a variety of adsorption properties for porous materials, including MOFs with open-metal sites and all-silica zeolites. Specifically, by using ESP-MMs, the binding geometry and adsorption energy landscape can be well captured. This enables these models to be employed to unravel the fundamental mechanism of gaseous adsorption in materials of interest as well as to facilitate the parametrization of adsorbent-adsorbate interactions. We also demonstrate that, combined with generic force fields for adsorbents, ESP-MMs can offer reasonable predictions in adsorption isotherms. Although these ESP-MMs use a relatively simple and nonpolarizable potential form for the sake of efficiency and applicability, their accuracy has been extensively validated in this study. Furthermore, the set of Lennard-Jones potentials with static point charges adopted for ESP-MMs can be readily implemented in all available simulation packages. We anticipate that these ESP-MMs can largely facilitate future computational studies of porous materials for gas separation and removal.

Ligand Binding Rate Constants in Heme Proteins Using Markov State Models and Molecular Dynamics Simulations

Computer simulation studies of the molecular basis for ligand migration in proteins allow the description of key events such as the transition between docking sites, displacement of existing ligands and solvent molecules, and open/closure of specific "gates", among others. In heme proteins, ligand migration from the solvent to the active site preludes the binding to the heme iron and triggers different functions. In this work, molecular dynamics simulations, a Markov State Model of migration and empirical kinetic equations are combined to study the migration of O₂ and NO in two truncated hemoglobins of Mycobacterium tuberculosis (Mt-TrHbN and Mt-TrHbO). For Mt-TrHbN, we show that the difference in the association constant in the oxy and deoxy states relies mainly in the displacement of water molecules anchored in the distal cavity in the deoxy form. The results here provide a valuable approach to study ligand migration in globins.

Large-Scale Computational Screening of Metal Organic Framework (MOF) Membranes and MOF-Based Polymer Membranes for H2/N2 Separations

Several thousands of metal organic frameworks (MOFs) have been reported to date, but the information on H₂/N₂ separation performances of MOF membranes is currently very limited in the literature. We report the first large-scale computational screening study that combines state-of-the-art molecular simulations, grand canonical Monte Carlo (GCMC) and molecular dynamics (MD), to predict H₂ permeability and H₂/N₂ selectivity of 3765 different types of MOF membranes. Results showed that MOF membranes offer very high H₂ permeabilities, 2.5 × 10³ to 1.7 × 10⁶ Barrer, and moderate H₂/N₂ membrane selectivities up to 7. The top 20 MOF membranes that exceed the polymeric membranes' upper bound for H₂/N₂ separation were identified based on the results of initial screening performed at infinite dilution condition. Molecular simulations were then carried out considering binary H₂/N₂ and quaternary H₂/N₂/CO₂/CO mixtures to evaluate the separation performance of MOF membranes under industrial operating conditions. Lower H₂ permeabilities and higher N₂ permeabilities were obtained at binary mixture conditions compared to the ones obtained at infinite dilution due to the absence of multicomponent mixture effects in the latter. Structure-performance relations of MOFs were also explored to provide molecular-level insights into the development of new MOF membranes that can offer both high H₂ permeability and high H₂/N₂ selectivity. Results showed that the most promising MOF membranes generally have large pore sizes (>6 Å) as well as high surface areas (>3500 m²/g) and high pore volumes (>1 cm³/g). We finally examined H₂/N₂ separation potentials of the mixed matrix membranes (MMMs) in which the best MOF materials identified from our high-throughput screening were used as fillers in various polymers. Results showed that incorporation of MOFs into polymers almost doubles H₂ permeabilities and slightly enhances H₂/N₂ selectivities of polymer membranes, which can advance the current membrane technology for efficient H₂ purification.

New Model for Predicting Adsorption of Polar Molecules in Metal-Organic Frameworks with Unsaturated Metal Sites

Conventional molecular models fail to correctly describe interactions of adsorbates with coordinatively unsaturated sites (CUS) present in a large number of metal-organic frameworks (MOFs). Here, we confirm the failure of these models for a prototypical polar adsorbate, carbon monoxide, and show that simply adjusting their parameters leads to poor agreement with experimental isotherms when outside the fitting conditions. We propose a new approach that combines quantum mechanical density functional theory (DFT) with Monte Carlo simulations to rigorously account for specific interactions at the CUS. By explicitly including electrostatic interactions and employing accurate DFT functionals that describe dispersion interactions, our modeling approach becomes generally applicable to both polar and nonpolar molecules. We demonstrate that this CUS model leads to substantial improvement in carbon monoxide adsorption isotherm predictions, and correctly captures the coordination binding mechanism. This paper represents a major stepping stone in the development of a robust, transferable and generally applicable approach to describe the complex interactions between gas molecules and CUS, with great potential for use in large-scale screening studies.

Kinetics of CO2 diffusion in human carbonic anhydrase: a study using molecular dynamics simulations and the Markov-state model

Molecular dynamics (MD) simulations, in combination with the Markov-state model (MSM), were applied to probe CO₂ diffusion from an aqueous solution into the active site of human carbonic anhydrase II (hCA-II), an enzyme useful for enhanced CO₂ capture and utilization. The diffusion process in the hydrophobic pocket of hCA-II was illustrated in terms of a two-dimensional free-energy landscape. We found that CO₂ diffusion in hCA-II is a rate-limiting step in the CO₂ diffusion-binding-reaction process. The equilibrium distribution of CO₂ shows its preferential accumulation within a hydrophobic domain in the protein core region. An analysis of the committors and reactive fluxes indicates that the main pathway for CO₂ diffusion into the active site of hCA-II is through a binding pocket where residue Gln¹³⁶ contributes to the maximal flux. The simulation results offer a new perspective on the CO₂ hydration kinetics and useful insights toward the development of novel biochemical processes for more efficient CO₂ sequestration and utilization.

Solvent Composition Drives the Rebinding Kinetics of Nitric Oxide to Microperoxidase

The rebinding kinetics of NO after photodissociation from microperoxidase (Mp-9) is studied in different solvent environments. In mixed glycerol/water (G/W) mixtures the dissociating ligand rebinds with a yield close to 1 due to the cavities formed by the solvent whereas in pure water the ligand can diffuse into the solvent after photodissociation. In the G/W mixture, only geminate rebinding on the sub-picosecond and 5 ps time scales was found and the rebinding fraction is unity which compares well with available experiments. Contrary to that, simulations in pure water find two time scales – ~10 ps and ~200 ps - indicating that both, geminate rebinding and rebinding after diffusion of NO in the surrounding water contribute. The rebinding fraction is around 0.63 within 1 ns which is in stark contrast with experiment. Including ions (Na and Cl) at 0.15 M concentration in water leads to rebinding kinetics tending to that in the glycerol/water mixture and yields agreement with experiments. The effect of temperature is also probed and found to be non-negligible. The present simulations suggest that NO rebinding in Mp is primarily driven by thermal fluctuations which is consistent with recent resonance Raman spectroscopy experiments and simulations on MbNO.

Ultrafast dynamics induced by the interaction of molecules with electromagnetic fields: Several quantum, semiclassical, and classical approaches

Several strategies for simulating the ultrafast dynamics of molecules induced by interactions with electromagnetic fields are presented. After a brief overview of the theory of molecule-field interaction, we present several representative examples of quantum, semiclassical, and classical approaches to describe the ultrafast molecular dynamics, including the multiconfiguration time-dependent Hartree method, Bohmian dynamics, local control theory, semiclassical thawed Gaussian approximation, phase averaging, dephasing representation, molecular mechanics with proton transfer, and multipolar force fields. In addition to the general overview, some focus is given to the description of nuclear quantum effects and to the direct dynamics, in which the ab initio energies and forces acting on the nuclei are evaluated on the fly. Several practical applications, performed within the framework of the Swiss National Center of Competence in Research "Molecular Ultrafast Science and Technology," are presented: These include Bohmian dynamics description of the collision of H with H2, local control theory applied to the photoinduced ultrafast intramolecular proton transfer, semiclassical evaluation of vibrationally resolved electronic absorption, emission, photoelectron, and time-resolved stimulated emission spectra, infrared spectroscopy of H-bonding systems, and multipolar force fields applications in the condensed phase.

A method for extracting carboxy-myoglobin from beef

The present study was conducted to devise a method for the effective extraction of carboxy-myoglobin (COMb) from beef without carbon monoxide dissociation. The ratio of COMb to myoglobin was computed at absorptions of wavelengths 541 and 551 nm, which characterize COMb and the isosbestic point between COMb and deoxy-myoglobin, respectively. The COMb extraction rate was found to vary with temperature, pH and oxygen conditions. The decrease observed in this rate was inversely proportional to the rise in extraction temperature. The COMb extraction rate was also affected by pH, and the stability of COMb in the extract solution was the highest at pH 8.0-9.0. Moreover, the presence of oxygen was found to disturb COMb extraction. According to these results, nearly all COMb could be extracted from carbon-monoxide-treated beef under stirring conditions in pH 8.5 deoxidized buffer, at 1°C, and under N₂ flow with the improved extraction method in this study (98.1 ± 2.7%). The decrement of COMb in the extract was accelerated by light, and the COMb was stable for 20 min in the dark, at 1°C. The extraction conditions for COMb described above should allow the accurate evaluation of COMb in meat tissue.

CO Separation from H2 via Hydrate Formation in Single-Walled Carbon Nanotubes

Hydrogen is an alternative fuel without generating greenhouse gas or other harmful emissions. Industrial hydrogen production, however, always contains a small fraction of carbon monoxide (CO) (∼0.5-2%) that must be removed for use in fuel cells. Here, we present molecular dynamics simulation evidence on facile separation of CO from H₂ at ambient pressure via the formation of quasi-one-dimensional (Q1D) clathrate hydrates within single-walled carbon nanotubes (SW-CNTs). At ambient pressure, Q1D CO (or H₂) clathrates in SW-CNTs are formed spontaneously when the SW-CNTs are immersed in CO (or H₂) aqueous solution. More interestingly, for the CO/H₂ aqueous solution, highly preferential adsorption of CO over H₂ occurs within the octagonal or nonagonal ice nanotubes inside of SW-CNTs. These results suggest that the formation of Q1D hydrates within SW-CNTs can be a viable and safe method for the separation of CO from H₂, which can be exploited for hydrogen purification in fuel cells.

The Dynamics Behind the Affinity: Controlling Heme-Gas Affinity via Geminate Recombination and Heme Propionate Conformation in the NO Carrier Cytochrome c'

Nitric oxide (NO) sensors are heme proteins which may also bind CO and O₂. Control of heme-gas affinity and their discrimination are achieved by the structural properties and reactivity of the heme and its distal and proximal environments, leading to several energy barriers. In the bacterial NO sensor cytochrome c' from Alcaligenes xylosoxidans (AXCP), the single Leu16Ala distal mutation boosts the affinity for gas ligands by a remarkable 10⁶-10⁸-fold, transforming AXCP from one of the lowest affinity gas binding proteins to one of the highest. Here, we report the dynamics of diatomics after photodissociation from wild type and L16A-AXCP over 12 orders of magnitude in time. For the L16A variant, the picosecond geminate rebinding of both CO and NO appears with an unprecedented 100% yield, and no exit of these ligands from protein to solvent could be observed. Molecular dynamic simulations saliently demonstrate that dissociated CO stays within 4 Å from Fe²⁺, in contrast to wild-type AXCP. The L16A mutation confers a heme propionate conformation and docking site which traps the diatomics, maximizing the probability of recombination and directly explaining the ultrahigh affinities for CO, NO, and O₂. Overall, our results point to a novel mechanism for modulating heme-gas affinities in proteins.

Ultrafast anisotropic protein quake propagation after CO photodissociation in myoglobin

"Protein quake" denotes the dissipation of excess energy across a protein, in response to a local perturbation such as the breaking of a chemical bond or the absorption of a photon. Femtosecond time-resolved small- and wide-angle X-ray scattering (TR-SWAXS) is capable of tracking such ultrafast protein dynamics. However, because the structural interpretation of the experiments is complicated, a molecular picture of protein quakes has remained elusive. In addition, new questions arose from recent TR-SWAXS data that were interpreted as underdamped oscillations of an entire protein, thus challenging the long-standing concept of overdamped global protein dynamics. Based on molecular-dynamics simulations, we present a detailed molecular movie of the protein quake after carbon monoxide (CO) photodissociation in myoglobin. The simulations suggest that the protein quake is characterized by a single pressure peak that propagates anisotropically within 500 fs across the protein and further into the solvent. By computing TR-SWAXS patterns from the simulations, we could interpret features in the reciprocal-space SWAXS signals as specific real-space dynamics, such as CO displacement and pressure wave propagation. Remarkably, we found that the small-angle data primarily detect modulations of the solvent density but not oscillations of the bare protein, thereby reconciling recent TR-SWAXS experiments with the notion of overdamped global protein dynamics.

Impact of Quadrupolar Electrostatics on Atoms Adjacent to the Sigma-Hole in Condensed-Phase Simulations

Halogenation is one of the cases for which advanced molecular simulation methods are mandatory for quantitative and predictive studies. The present work provides a systematic investigation of the importance of higher-order multipoles on specific sites of halobenzenes, other than the halogen, for static and dynamic properties in condensed-phase simulations. For that purpose, solute-solvent interactions using point charge (PC), multipole (MTP), and hybrid point charge/multipole (HYB) electrostatic models are analyzed in regions of halogen bonding and extended to regions of π orbitals of phenyl carbons. Using molecular dynamics simulations and quantum chemical methods, it is found that the sigma-hole does not only affect the halogen and the carbon bound to it but its effect extends to the carbons adjacent to the CX bond. This effect increases with the magnitude of the positive potential of the sigma-hole. With the MTP and HYB3 models, all hydration free energies of the PhX compounds are reproduced within 0.1 kcal/mol. Analysis of pair distribution functions and hydration free energies of halogenated benzenes provides a microscopic explanation why "point charge"-based representations with off-site charges fail in reproducing thermodynamic properties of the sigma-hole. Application of the hybrid models to study protein-ligand binding demonstrates both their accuracy and computational efficiency.

Nanosecond ligand migration and functional protein relaxation in ba3 oxidoreductase: Structures of the B0, B1 and B2 intermediate states

Nanosecond time-resolved step-scan FTIR spectroscopy (nTRS (2) -FTIR) has been applied to literally probe the active site of the carbon monoxide (CO)-bound thermophilic ba3 heme-copper oxidoreductase as it executes its function. The nTRS (2) - snapshots of the photolysed heme a3 Fe-CO/CuB species captured a "transition state" whose side chains prevent the photolysed CO to enter the docking cavity. There are three sets of ba3 photoproduct bands of docked CO with different orientation exhibiting different kinetics. The trajectories of the "docked" CO at 2122, 2129 and 2137cm(-1) is referred to in the literature as B2, B1 and B0 intermediate states, respectively. The present data provided direct evidence for the role of water in controlling ligand orientation in an intracavity protein environment.

Atomistic Simulation of Adiabatic Reactive Processes Based on Multi-State Potential Energy Surfaces

The adiabatic reactive molecular dynamics (ARMD) method provides a framework to study chemical reactions using molecular dynamics simulations with minimal computational overhead. Here, ARMD is generalized to an arbitrary reactive process between two states in which reactants and products can be treated by an atomistic force field. The implementation is described, and the method is applied to two systems: the kinetics of NO rebinding to myoglobin (Mb) as a validation system and the conformational transition in neuroglobin (Ngb) which explores the full functionality of ARMD. For MbNO, the nonexponential kinetics observed both in experiment and earlier ARMD studies is reproduced. Furthermore, the sensitivity of the results with respect to the asymptotic separation between the two potential energy surfaces (NO bound and unbound) is studied.

Identification of Mutational Hot Spots for Substrate Diffusion: Application to Myoglobin

The pathways by which small molecules (substrates or inhibitors) access active sites are a key aspect of the function of enzymes and other proteins. A key problem in designing or altering such proteins is to identify sites for mutation that will have the desired effect on the substrate transport properties. While specific access channels have been invoked in the past, molecular simulations suggest that multiple routes are possible, complicating the analysis. This complexity, however, can be captured by a Markov State Model (MSM) of the ligand diffusion process. We have developed a sensitivity analysis of the resulting rate matrix, which identifies the locations where mutations should have the largest effect on the diffusive on rate. We apply this method to myoglobin, which is the best characterized example both from experiment and simulation. We validate the approach by translating the sensitivity parameter obtained from this method into the CO binding rates in myoglobin upon mutation, resulting in a semi-quantitative correlation with experiments. The model is further validated against an explicit simulation for one of the experimental mutants.

Diffusion of molecules in the bulk of a low density amorphous ice from molecular dynamics simulations

Arguments for a solvent driven mechanism for the diffusion of CO, CO₂, NH₃, and H₂CO in a LDA water ice.

Non-site-specific allosteric effect of oxygen on human hemoglobin under high oxygen partial pressure

Protein allostery is essential for vital activities. Allosteric regulation of human hemoglobin (HbA) with two quaternary states T and R has been a paradigm of allosteric structural regulation of proteins. It is widely accepted that oxygen molecules (O₂) act as a “site-specific” homotropic effector, or the successive O₂ binding to the heme brings about the quaternary regulation. However, here we show that the site-specific allosteric effect is not necessarily only a unique mechanism of O₂ allostery. Our simulation results revealed that the solution environment of high O₂ partial pressure enhances the quaternary change from T to R without binding to the heme, suggesting an additional “non-site-specific” allosteric effect of O₂. The latter effect should play a complementary role in the quaternary change by affecting the intersubunit contacts. This analysis must become a milestone in comprehensive understanding of the allosteric regulation of HbA from the molecular point of view.

2D IR spectra of cyanide in water investigated by molecular dynamics simulations

Using classical molecular dynamics simulations, the 2D infrared (IR) spectroscopy of CN(-) solvated in D2O is investigated. Depending on the force field parametrizations, most of which are based on multipolar interactions for the CN(-) molecule, the frequency-frequency correlation function and observables computed from it differ. Most notably, models based on multipoles for CN(-) and TIP3P for water yield quantitatively correct results when compared with experiments. Furthermore, the recent finding that T1 times are sensitive to the van der Waals ranges on the CN(-) is confirmed in the present study. For the linear IR spectrum, the best model reproduces the full widths at half maximum almost quantitatively (13.0 cm(-1) vs. 14.9 cm(-1)) if the rotational contribution to the linewidth is included. Without the rotational contribution, the lines are too narrow by about a factor of two, which agrees with Raman and IR experiments. The computed and experimental tilt angles (or nodal slopes) α as a function of the 2D IR waiting time compare favorably with the measured ones and the frequency fluctuation correlation function is invariably found to contain three time scales: a sub-ps, 1 ps, and one on the 10-ps time scale. These time scales are discussed in terms of the structural dynamics of the surrounding solvent and it is found that the longest time scale (≈10 ps) most likely corresponds to solvent exchange between the first and second solvation shell, in agreement with interpretations from nuclear magnetic resonance measurements.

Oxygen migration pathways in NO-bound truncated hemoglobin

Atomistic simulations of dioxygen (O(2)) dynamics and migration in nitric oxide-bound truncated Hemoglobin N (trHbN) of Mycobacterium tuberculosis are reported. From more than 100 ns of simulations the connectivity network involving the metastable states for localization of the O(2) ligand is built and analyzed. It is found that channel I is the primary entrance point for O(2) whereas channel II is predominantly an exit path although access to the protein active site is also possible. For O(2) a new site compared to nitric oxide, from which reaction with the heme group can occur, was found. As this site is close to the heme iron, it could play an important role in the dioxygenation mechanism as O(2) can remain there for hundreds of picoseconds after which it can eventually leave the protein, while NO is localized in Xe2. The present study supports recent experimental work which proposed that O(2) docks in alternative pockets than Xe close to the reactive site. Similar to other proteins, a phenylalanine residue (Phe62) plays the role of a gate along the access route in channel I. The most highly connected site is the Xe3 pocket which is a "hub" and free energy barriers between the different metastable states are ≈1.5 kcal mol(-1) which allows facile O(2) migration within the protein.

Ligand migration in myoglobin: a combined study of computer simulation and x-ray crystallography

A ligand-migration mechanism of myoglobin was studied by a multidisciplinary approach that used x-ray crystallography and molecular dynamics simulation. The former revealed the structural changes of the protein along with the ligand migration, and the latter provided the statistical ensemble of protein conformations around the thermal average. We developed a novel computational method, homogeneous ensemble displacement, and generated the conformational ensemble of ligand-detached species from that of ligand-bound species. The thermally averaged ligand-protein interaction was illustrated in terms of the potential of mean force. Although the structural changes were small, the presence of the ligand molecule in the protein matrix significantly affected the 3D scalar field of the potential of mean force, in accordance with the self-opening model proposed in the previous x-ray study.

Quantifying the importance of protein conformation on ligand migration in myoglobin

Myoglobin (Mb) is a model system for ligand binding and migration. The energy barriers (ΔG) for ligand migration in Mb have been studied in the past by experiment and theory and significant differences between different approaches were found. From experiment, it is known that Mb can assume a large number of conformational substates. In this work, these substates are investigated as a possible source of the differences in migration barriers. We show that the initial structure significantly affects the calculated ΔG for a particular transition and that fluctuations in barrier heights δΔG are of similar magnitude as the free energy barriers themselves. The sensitivity of ΔG to the initial structure is compared to other sources of errors. Different protein structures can affect the calculated ΔG by up to 4 kcal/mol, whereas differences between simple point charge models and more elaborate multipolar charge models for the CO-ligand are smaller by a factor of two. Analysis of the structural changes underlying the large effect of the conformational substate reveals the importance of coupling between protein and ligand motion for migration.

Kinetics of carbon monoxide migration and binding in solvated neuroglobin as revealed by molecular dynamics simulations and quantum mechanical calculations

Neuroglobin (Ngb) is a globular protein that reversibly binds small ligands at the six coordination position of the heme. With respect to other globins similar to myoglobin, Ngb displays some peculiarities as the topological reorganization of the internal cavities coupled to the sliding of the heme, or the binding of the endogenous distal histidine to the heme in the absence of an exogenous ligand. In this Article, by using multiple (independent) molecular dynamics trajectories (about 500 ns in total), the migration pathways of photolized carbon monoxide (CO) within solvated Ngb were analyzed, and a quantitative description of CO migration and corresponding kinetics was obtained. MD results, combined with quantum mechanical calculations on the CO-heme binding-unbinding reaction step in Ngb, allowed construction of a quantitative model representing the relevant steps of CO migration and rebinding.

Design of an infrared laser pulse to control the multiphoton dissociation of the Fe-CO bond in CO-heme compounds

Optimal control theory is used to design a laser pulse for the multiphoton dissociation of the Fe-CO bond in the CO-heme compounds. The study uses a hexacoordinated iron-porphyrin-imidazole-CO complex in its ground electronic state as a model for CO liganded to the heme group. The potential energy and dipole moment surfaces for the interaction of the CO ligand with the heme group are calculated using density functional theory. Optimal control theory, combined with a time-dependent quantum dynamical treatment of the laser-molecule interaction, is then used to design a laser pulse capable of efficiently dissociating the CO-heme complex model. The genetic algorithm method is used within the mathematical framework of optimal control theory to perform the optimization process. This method provides good control over the parameters of the laser pulse, allowing optimized pulses with simple time and frequency structures to be designed. The dependence of photodissociation yield on the choice of initial vibrational state and of initial laser field parameters is also investigated. The current work uses a reduced dimensionality model in which only the Fe-C and C-O stretching coordinates are explicitly taken into account in the time-dependent quantum dynamical calculations. The limitations arising from this are discussed in Sec. IV.