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

SCN as a local probe of protein structural dynamics

The dynamics of lysozyme is probed by attaching -SCN to all alanine residues. The one-dimensional infrared spectra exhibit frequency shifts in the position of the maximum absorption of 4 cm-1, which is consistent with experiments in different solvents and indicates moderately strong interactions of the vibrational probe with its environment. Isotopic substitution 12C → 13C leads to a redshift by -47 cm-1, which agrees quantitatively with experiments for CN-substituted copper complexes in solution. The low-frequency, far-infrared part of the protein spectra contains label-specific information in the difference spectra when compared with the wild type protein. Depending on the position of the labels, local structural changes are observed. For example, introducing the -SCN label at Ala129 leads to breaking of the α-helical structure with concomitant change in the far-infrared spectrum. Finally, changes in the local hydration of SCN-labeled alanine residues as a function of time can be related to the reorientation of the label. It is concluded that -SCN is potentially useful for probing protein dynamics, both in the high-frequency part (CN-stretch) and in the far-infrared part of the spectrum.

Development of a universal method for vibrational analysis of the terminal alkyne C≡C stretch

The terminal alkyne C≡C stretch has a large Raman scattering cross section in the "silent" region for biomolecules. This has led to many Raman tag and probe studies using this moiety to study biomolecular systems. A computational investigation of these systems is vital to aid in the interpretation of these results. In this work, we develop a method for computing terminal alkyne vibrational frequencies and isotropic transition polarizabilities that can easily and accurately be applied to any terminal alkyne molecule. We apply the discrete variable representation method to a localized version of the C≡C stretch normal mode. The errors of (1) vibrational localization to the terminal alkyne moiety, (2) anharmonic normal mode isolation, and (3) discretization of the Born-Oppenheimer potential energy surface are quantified and found to be generally small and cancel each other. This results in a method with low error compared to other anharmonic vibrational methods like second-order vibrational perturbation theory and to experiments. Several density functionals are tested using the method, and TPSS-D3, an inexpensive nonempirical density functional with dispersion corrections, is found to perform surprisingly well. Diffuse basis functions are found to be important for the accuracy of computed frequencies. Finally, the computation of vibrational properties like isotropic transition polarizabilities and the universality of the localized normal mode for terminal alkynes are demonstrated.

Structural basis for bacterial energy extraction from atmospheric hydrogen

Diverse aerobic bacteria use atmospheric H₂ as an energy source for growth and survival¹. This globally significant process regulates the composition of the atmosphere, enhances soil biodiversity and drives primary production in extreme environments^2,3. Atmospheric H₂ oxidation is attributed to uncharacterized members of the [NiFe] hydrogenase superfamily^4,5. However, it remains unresolved how these enzymes overcome the extraordinary catalytic challenge of oxidizing picomolar levels of H₂ amid ambient levels of the catalytic poison O₂ and how the derived electrons are transferred to the respiratory chain¹. Here we determined the cryo-electron microscopy structure of the Mycobacterium smegmatis hydrogenase Huc and investigated its mechanism. Huc is a highly efficient oxygen-insensitive enzyme that couples oxidation of atmospheric H₂ to the hydrogenation of the respiratory electron carrier menaquinone. Huc uses narrow hydrophobic gas channels to selectively bind atmospheric H₂ at the expense of O₂, and 3 [3Fe-4S] clusters modulate the properties of the enzyme so that atmospheric H₂ oxidation is energetically feasible. The Huc catalytic subunits form an octameric 833 kDa complex around a membrane-associated stalk, which transports and reduces menaquinone 94 Å from the membrane. These findings provide a mechanistic basis for the biogeochemically and ecologically important process of atmospheric H₂ oxidation, uncover a mode of energy coupling dependent on long-range quinone transport, and pave the way for the development of catalysts that oxidize H₂ in ambient air.

Structure, Organization, and Heterogeneity of Water-Containing Deep Eutectic Solvents

The spectroscopy and structural dynamics of a deep eutectic mixture (KSCN/acetamide) with varying water content is investigated from 2D IR (with the C-N stretch vibration of the SCN^- anions as the reporter) and THz spectroscopy. Molecular dynamics simulations correctly describe the nontrivial dependence of both spectroscopic signatures depending on water content. For the 2D IR spectra, the MD simulations relate the steep increase in the cross-relaxation rate at high water content to the parallel alignment of packed SCN^- anions. Conversely, the nonlinear increase of the THz absorption with increasing water content is mainly attributed to the formation of larger water clusters. The results demonstrate that a combination of structure-sensitive spectroscopies and molecular dynamics simulations provides molecular-level insights into the emergence of heterogeneity of such mixtures by modulating their composition.

Quantitative molecular simulations

All-atom simulations can provide molecular-level insights into the dynamics of gas-phase, condensed-phase and surface processes. One important requirement is a sufficiently realistic and detailed description of the underlying intermolecular interactions. The present perspective provides an overview of the present status of quantitative atomistic simulations from colleagues' and our own efforts for gas- and solution-phase processes and for the dynamics on surfaces. Particular attention is paid to direct comparison with experiment. An outlook discusses present challenges and future extensions to bring such dynamics simulations even closer to reality.

Investigating the effects of N-terminal acetylation on KFE8 self-assembly with 2D IR spectroscopy

Peptide self-assembly is an exciting and robust approach to create novel nanoscale materials for biomedical applications. However, the complex interplay between intra- and intermolecular interactions in peptide aggregation means that minor changes in peptide sequence can yield dramatic changes in supramolecular structure. Here, we use two-dimensional infrared (2D IR) spectroscopy to study a model amphiphilic peptide, KFE8, and its N-terminal acetylated counterpart, AcKFE8. 2D IR spectra of isotope-labeled peptides reveal that AcKFE8 aggregates comprise two distinct β-sheet structures while KFE8 aggregates comprise only one of these structures. Using an excitonic Hamiltonian to simulate the vibrational spectra of model β-sheets, we determine that the spectra are consistent with antiparallel β-sheets with different strand alignments, specifically a two-residue shift in the register of the β-strands. These findings bring forth new insights into how N-terminal acetylation may subtly impact secondary structure, leading to larger effects on overall aggregate morphology. Additionally, these results highlight the importance of understanding the residue-level structural differences that result from changes in peptide sequence in order to facilitate the rational design of peptide materials.

Hydration dynamics and IR spectroscopy of 4-fluorophenol

Halogenated groups are relevant in pharmaceutical applications and potentially useful spectroscopic probes for infrared spectroscopy. In this work, the structural dynamics and infrared spectroscopy of para-fluorophenol (F-PhOH) and phenol (PhOH) is investigated in the gas phase and in water using a combination of experiment and molecular dynamics (MD) simulations. The gas phase and solvent dynamics around F-PhOH and PhOH is characterized from atomistic simulations using empirical energy functions with point charges or multipoles for the electrostatics, Machine Learning (ML) based parametrizations and with full ab initio (QM) and mixed Quantum Mechanical/Molecular Mechanics (QM/MM) simulations with a particular focus on the CF- and OH-stretch region. The CF-stretch band is heavily mixed with other modes whereas the OH-stretch in solution displays a characteristic high-frequency peak around 3600 cm⁻¹ most likely associated with the –OH group of PhOH and F-PhOH together with a characteristic progression below 3000 cm⁻¹ due to coupling with water modes which is also reproduced by several of the simulations. Solvent and radial distribution functions indicate that the CF-site is largely hydrophobic except for simulations using point charges which renders them unsuited for correctly describing hydration and dynamics around fluorinated sites. The hydrophobic character of the CF-group is particularly relevant for applications in pharmaceutical chemistry with a focus on local hydration and interaction with the surrounding protein.

Halogenated groups are relevant in pharmaceutical applications and potentially useful spectroscopic probes for infrared spectroscopy.

Multipolar Force Fields for Amide-I Spectroscopy from Conformational Dynamics of the Alanine Trimer

The dynamics and spectroscopy of N-methyl-acetamide (NMA) and trialanine in solution are characterized from molecular dynamics simulations using different energy functions, including a conventional point charge (PC)-based force field, one based on a multipolar (MTP) representation of the electrostatics, and a semiempirical DFT method. For the 1D infrared spectra, the frequency splitting between the two amide-I groups is 10 cm^-1 from the PC, 13 cm^-1 from the MTP, and 47 cm^-1 from self-consistent charge density functional tight-binding (SCC-DFTB) simulations, compared with 25 cm^-1 from experiment. The frequency trajectory required for the frequency fluctuation correlation function (FFCF) is determined from individual normal mode (INM) and full normal mode (FNM) analyses of the amide-I vibrations. The spectroscopy, time-zero magnitude of the FFCF C(t = 0), and the static component Δ₀² from simulations using MTP and analysis based on FNM are all consistent with experiments for (Ala)₃. Contrary to this, for the analysis excluding mode-mode coupling (INM), the FFCF decays to zero too rapidly and for simulations with a PC-based force field, the Δ₀² is too small by a factor of two compared with experiments. Simulations with SCC-DFTB agree better with experiment for these observables than those from PC-based simulations. The conformational ensemble sampled from simulations using PCs is consistent with the literature (including P_II, β, α_R, and α_L), whereas that covered by the MTP-based simulations is dominated by P_II with some contributions from β and α_R. This agrees with and confirms recently reported Bayesian-refined populations based on 1D infrared experiments. FNM analysis together with a MTP representation provides a meaningful model to correctly describe the dynamics of hydrated trialanine.

Site-selective dynamics of azidolysozyme

The spectroscopic response of and structural dynamics around all azido-modified alanine residues (AlaN₃) in lysozyme are characterized. It is found that AlaN₃ is a positionally sensitive probe for the local dynamics, covering a frequency range of ∼15 cm^-1 for the center frequency of the line shape. This is consistent with findings from selective replacements of amino acids in PDZ2, which reported a frequency span of ∼10 cm^-1 for replacements of Val, Ala, or Glu by azidohomoalanine. For the frequency fluctuation correlation functions, the long-time decay constants τ₂ range from ∼1 to ∼10 ps, which compares with experimentally measured correlation times of 3 ps. Attaching azide to alanine residues can yield dynamics that decays to zero on the few ps time scale (i.e., static component Δ₀ ∼ 0 ps^-1) or to a remaining, static contribution of ∼0.5 ps^-1 (corresponding to 2.5 cm^-1), depending on the local environment on the 10 ps time scale. The magnitude of the static component correlates qualitatively with the degree of hydration of the spectroscopic probe. Although attaching azide to alanine residues is found to be structurally minimally invasive with respect to the overall protein structure, analysis of the local hydrophobicity indicates that the hydration around the modification site differs for modified and unmodified alanine residues, respectively.

Dynamics and Infrared Spectrocopy of Monomeric and Dimeric Wild Type and Mutant Insulin

The infrared spectroscopy and dynamics of -CO labels in wild type and mutant insulin monomer and dimer are characterized from molecular dynamics simulations using validated force fields. It is found that the spectroscopy of monomeric and dimeric forms in the region of the amide-I vibration differs for residues B24-B26 and D24-D26, which are involved in dimerization of the hormone. Also, the spectroscopic signatures change for mutations at position B24 from phenylalanine, which is conserved in many organisms and is known to play a central role in insulin aggregation, to alanine or glycine. Using three different methods to determine the frequency trajectories (solving the nuclear Schrödinger equation on an effective 1-dimensional potential energy curve, using instantaneous normal modes, and using parametrized frequency maps) leads to the same overall conclusions. The spectroscopic response of monomeric WT and mutant insulin differs from that of their respective dimers, and the spectroscopy of the two monomers in the dimer is also not identical. For the WT and F24A and F24G monomers, spectroscopic shifts are found to be ∼20 cm^-1 for residues (B24-B26) located at the dimerization interface. Although the crystal structure of the dimer is that of a symmetric homodimer, dynamically the two monomers are not equivalent on the nanosecond time scale. Together with earlier work on the thermodynamic stability of the WT and the same mutants, it is concluded that combining computational and experimental infrared spectroscopy provides a potentially powerful way to characterize the aggregation state and dimerization energy of modified insulins.

Multistate Reactive Molecular Dynamics Simulations of Proton Diffusion in Water Clusters and in the Bulk

The molecular mechanics with proton transfer (MMPT) force field is combined with multistate adiabatic reactive molecular dynamics (MS-ARMD) to describe proton transport in the condensed phase. Parametrization for small protonated water clusters based on electronic structure calculations at the MP2/6-311+G(2d,2p) level of theory and refinement by comparing with infrared spectra for a protonated water tetramer yields a force field which faithfully describes the minimum energy structures of small protonated water clusters. In protonated water clusters up to (H₂O)₁₀₀H⁺, the proton hopping rate is around 100 hops/ns. This rate converges for 21 ≤ n ≤ 31, and no further speedup in bulk water is found. This indicates that bulklike behavior requires the solvation of a Zundel motif by ∼25 water molecules, which corresponds to the second solvation sphere. For smaller cluster sizes, the number of available states (i.e., the number of proton acceptors) is too small and slows down proton-transfer rates. The cluster simulations confirm that the excess proton is typically located on the surface. The free-energy surface as a function of the weights of the two lowest states and a configurational parameter suggests that the "special pair" plays a central role in rapid proton transport. The barriers between this minimum-energy structure and the Zundel and Eigen minima are sufficiently low (∼1 kcal/mol, consistent with recent experiments and commensurate with a hopping rate of ∼100/ns or 1 every 10 ps), leading to a highly dynamic environment. These findings are also consistent with recent experiments which find that Zundel-type hydration geometries are prevalent in bulk water.

Dynamics of linear molecules in water: Translation-rotation coupling in jump motion driven diffusion

We study by computer simulations, and by theory, the coupled rotational and translational dynamics of three important linear diatomic molecules, namely, carbon monoxide (CO), nitric oxide (NO), and cyanide ion (CN^-) in water. Translational diffusion of these molecules is found to be strongly coupled to their own rotational dynamics which, in turn, are coupled to similar motions of the surrounding water. In particular, we find that coupled orientational jump motions play an important role in all three cases. While CO and NO show similar features, CN^- exhibits certain differences. Our results agree well with the known experimental values of the diffusion coefficient. We examined the validity of hydrodynamic predictions and found them to be inadequate, particularly for rotational diffusion. A mode coupling theory approach is developed and applied to understand the complexity of translation-rotation coupling.

Probing the Differential Dynamics of the Monomeric and Dimeric Insulin from Amide-I IR Spectroscopy

The monomer-dimer equilibrium for insulin is one of the essential steps in forming the receptor-binding competent monomeric form of the hormone. Despite this importance, the thermodynamic stability, in particular for modified insulins, is quite poorly understood, in part, due to experimental difficulties. This work explores one- and two-dimensional infrared spectroscopy in the range of the amide-I band for the hydrated monomeric and dimeric wild-type hormone. It is found that for the monomer the frequency fluctuation correlation function (FFCF) and the one-dimensional infrared spectra are position sensitive. The spectra of the -CO probes at the dimerization interface (residues Phe24, Phe25, and Tyr26) split and indicate an asymmetry despite the overall (formal) point symmetry of the dimer structure. Also, the decay times of the FFCF for the same -CO probe in the monomer and the dimer can differ by up to 1 order of magnitude, for example, for residue PheB24, which is solvent exposed for the monomer but at the interface for the dimer. The spectroscopic shifts correlate approximately with the average number of hydration waters and the magnitude of the FFCF at time zero. However, this correlation is only qualitative due to the heterogeneous and highly dynamical environment. Based on density functional theory calculations, the dominant contribution for solvent-exposed -CO is found to arise from the surrounding water (∼75%), whereas the protein environment contributes considerably less. The results suggest that infrared spectroscopy is a positionally sensitive probe of insulin dimerization, in particular in conjunction with isotopic labeling of the probe.

Vibrational Spectroscopy of N3- in the Gas and Condensed Phase

Azido-derivatized amino acids are potentially useful, positionally resolved spectroscopic probes for studying the structural dynamics of proteins and macromolecules in solution. To this end, a computational model for the vibrational modes of N₃^- based on accurate electronic structure calculations and a reproducing kernel Hilbert space representation of the potential energy surface for the internal degrees of freedom is developed. Fully dimensional quantum bound state calculations yield the antisymmetric stretch vibration at 1974 cm^-1 compared with 1986 cm^-1 from experiment. This mode shifts by 64 cm^-1 (from the frequency distribution) and 74 cm^-1 (from the IR line shape) to the blue, respectively, compared with 61 cm^-1 from experiment for N₃^- in water. The decay time of the frequency fluctuation correlation function is 1.1 ps, which is in good agreement with experiment (1.2-1.3 ps) and the full width at half maximum of the asymmetric stretch in solution is 18.5 cm^-1 compared with 25.2 cm^-1 from experiment. A computationally more efficient analysis based on instantaneous normal modes is shown to provide comparable, albeit somewhat less quantitative results compared to solving the three-dimensional Schrödinger equation for the fundamental vibrations.

Unmasking Rare, Large-Amplitude Motions in D2-Tagged I-·(H2O)2 Isotopomers with Two-Color, Infrared-Infrared Vibrational Predissociation Spectroscopy

We describe a two-color, isotopomer-selective infrared-infrared population-labeling method that can monitor very slow spectral diffusion of OH oscillators in H-bonded networks and apply it to the I^-·(HDO)·(D₂O) and I^-·(H₂O)·(D₂O) systems, which are cryogenically cooled and D₂-tagged at an ion trap temperature of 15 K. These measurements reveal very large (>400 cm^-1), spontaneous spectral shifts despite the fact that the predissociation spectra in the OH stretching region of both isotopologues are sharp and readily assigned to four fundamentals of largely decoupled OH oscillators held in a cyclic H-bonded network. This spectral diffusion is not observed in the untagged isotopologues of the dihydrate clusters that are generated under the same source conditions but does become apparent at about 75 K. These results are discussed in the context of the large-amplitude "jump" mechanism for H-bond relaxation dynamics advanced by Laage and Hynes in an experimental scenario where rare events can be captured by following the migration of OH groups among the four available positions in the quasi-rigid equilibrium structure.

Minimal distributed charges: Multipolar quality at the cost of point charge electrostatics

Most empirical force fields use atom-centered point charges (PCs) to represent the electrostatic potential (ESP) around molecules. While such PC models are computationally efficient, they are unable to capture anisotropic electronic features, such as σ holes or lone pairs. These features are better described using atomic multipole (MTP) moments, which significantly improve the quality of the resulting ESP. However, the improvement comes at the expense of a considerably increased computational complexity and cost for calculating the interaction energies and forces. In the present work, a novel minimal distributed charge model (MDCM) based on off-centered point charges is presented and the quality of the resulting ESP is compared to the performance of MTPs and atom-centered PC models for several test molecules. All three models are fitted using the same algorithm based on differential evolution, which is available as a Fortran90 program from the authors upon request. We show that the MDCM is capable of approximating the reference ab initio ESP with an accuracy as good as, or better than, MTPs without the need for computationally expensive higher order multipoles. Further it is demonstrated that the MDCM is numerically stable in molecular dynamics simulations and is able to reproduce electrostatic interaction energies and thermodynamic quantities with the same accuracy as MTPs at reduced computational cost.

Vibrational Stark spectroscopy for assessing ligand-binding strengths in a protein

Nitrile groups are potentially useful spectroscopic probes in the infrared to characterize the binding and dynamics of ligands in proteins. This opens the possibility of locating and determining the binding mode of suitably labelled ligands in proteins based on optical spectroscopy, without the need for determining an X-ray structure. However, relating structure and spectroscopy requires means to accurately compute infrared spectra. This is investigated for benzonitrile (PhCN) in water, wild type (WT) and two lysozyme mutants in solution. The force field is validated by comparing with experimental data for benzonitrile in water which is the basis for computing the Stark shift and time scale for spectral diffusion of PhCN in WT and the L99A and L99G mutants of T4 lysozyme. The 1-d spectra for PhCN in WT and the two mutant proteins differ in their maximum absorption by up to 4 cm^-1, which reflects the modified electrostatic environments in the three proteins. It is also tested whether extending from 1-d to 2-d infrared spectroscopy provides further discrimination in the ligand-binding modes. First, for PhCN in solution the frequency fluctuation correlation function (FFCF) decays to zero at short times whereas in the protein a pronounced static inhomogeneous component is found. Secondly, the decay time of the FFCF for the mutant to which PhCN binds most strongly has the longest decay time. It is demonstrated explicitly that the ligand-binding free energy with respect to the three protein variants correlates with the Stark shift. This makes 1-d infrared spectroscopy together with computations a valuable tool for characterizing binding modes and potentially binding locations in proteins.

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.

Structure of Penta-Alanine Investigated by Two-Dimensional Infrared Spectroscopy and Molecular Dynamics Simulation

We have studied the structure of (Ala)5, a model unfolded peptide, using a combination of 2D IR spectroscopy and molecular dynamics (MD) simulation. Two different isotopomers, each bis-labeled with (13)C═O and (13)C═(18)O, were strategically designed to shift individual site frequencies and uncouple neighboring amide-I' modes. 2D IR spectra taken under the double-crossed ⟨π/4, -π/4, Y, Z⟩ polarization show that the labeled four-oscillator systems can be approximated by three two-oscillator systems. By utilizing the different polarization dependence of diagonal and cross peaks, we extracted the coupling constants and angles between three pairs of amide-I' transition dipoles through spectral fitting. These parameters were related to the peptide backbone dihedral angles through DFT calculated maps. The derived dihedral angles are all located in the polyproline-II (ppII) region of the Ramachandran plot. These results were compared to the conformations sampled by Hamiltonian replica-exchange MD simulations with three different CHARMM force fields. The C36 force field predicted that ppII is the dominant conformation, consistent with the experimental findings, whereas C22/CMAP predicted similar population for α+, β, and ppII, and the polarizable Drude-2013 predicted dominating β structure. Spectral simulation based on MD representative conformations and structure ensembles demonstrated the need to include multiple 2D spectral features, especially the cross-peak intensity ratio and shape, in structure determination. Using 2D reference spectra defined by the C36 structure ensemble, the best spectral simulation is achieved with nearly 100% ppII population, although the agreement with the experimental cross-peak intensity ratio is still insufficient. The dependence of population determination on the choice of reference structures/spectra and the current limitations on theoretical modeling relating peptide structures to spectral parameters are discussed. Compared with the previous results on alanine based oligopeptides, the dihedral angles of our fitted structure, and the most populated ppII structure from the C36 simulation are in good agreement with those suggesting a major ppII population. Our results provide further support for the importance of ppII conformation in the ensemble of unfolded peptides.

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.

Transferable Atomic Multipole Machine Learning Models for Small Organic Molecules

Accurate representation of the molecular electrostatic potential, which is often expanded in distributed multipole moments, is crucial for an efficient evaluation of intermolecular interactions. Here we introduce a machine learning model for multipole coefficients of atom types H, C, O, N, S, F, and Cl in any molecular conformation. The model is trained on quantum-chemical results for atoms in varying chemical environments drawn from thousands of organic molecules. Multipoles in systems with neutral, cationic, and anionic molecular charge states are treated with individual models. The models' predictive accuracy and applicability are illustrated by evaluating intermolecular interaction energies of nearly 1,000 dimers and the cohesive energy of the benzene crystal.

Ultrafast Dynamics of Carboxy-Hemoglobin: Two-Dimensional Infrared Spectroscopy Experiments and Simulations

This Letter presents a comparison between experimental and simulated 2D mid-infrared spectra of carboxy-hemoglobin in the spectral region of the carbon monoxide stretching mode. The simulations rely on a fluctuating potential energy surface that includes both the effect of heme and the protein surroundings computed from molecular dynamics simulations. A very good agreement between theory and experiment is obtained with no adjustable parameters. The simulations show that the effect of the distal histidine through the hydrogen bond is strong and is directly responsible for the slow decay of the frequency-frequency correlation function on a 10 ps time scale. This study confirms that fluctuations in carboxy-hemoglobin are more inhomogeneous than those in the more frequently studied carboxy-myoglobin. The comparison between simulations and experiments brings valuable information on the complex relation between protein structure and spectral diffusion.

Multipolar force fields and their effects on solvent dynamics around simple solutes

The performance of multipole (MTP) and point charge (PC) force fields in classical molecular dynamics (MD) simulations of condensed-phase systems for both equilibrium and dynamical quantities is compared. MTP electrostatics provides an improved description of the anisotropic electrostatic potential, which is especially important to describe key, challenging interactions, such as lone pairs, π-interactions, and hydrogen bonds. These chemical environments are probed by focusing on the hydration properties of two molecules: N-methylacetamide and phenyl bromide. Both, equilibrium and dynamical, quantities are affected by the quality of the electrostatic model. The alteration of the first solvation shell in MTP simulations is validated by comparing with lifetimes and correlation times of solute-solvent interactions from experiment. The improved dynamical behavior found in the MTP simulations—observed for molecules parametrized using very different protocols—suggests that a systematic improvement of both equilibrium and dynamical quantities when using MTP electrostatics is possible.

A Novel, Computationally Efficient Multipolar Model Employing Distributed Charges for Molecular Dynamics Simulations

A truncated multipole expansion can be re-expressed exactly using an appropriate arrangement of point charges. This means that groups of point charges that are shifted away from nuclear coordinates can be used to achieve accurate electrostatics for molecular systems. We introduce a multipolar electrostatic model formulated in this way for use in computationally efficient multipolar molecular dynamics simulations with well-defined forces and energy conservation in NVE (constant number-volume-energy) simulations. A framework is introduced to distribute torques arising from multipole moments throughout a molecule, and a refined fitting approach is suggested to obtain atomic multipole moments that are optimized for accuracy and numerical stability in a force field context. The formulation of the charge model is outlined as it has been implemented into CHARMM, with application to test systems involving H2O and chlorobenzene. As well as ease of implementation and computational efficiency, the approach can be used to provide snapshots for multipolar QM/MM calculations in QM/MM-MD studies and easily combined with a standard point-charge force field to allow mixed multipolar/point charge simulations of large systems.

Coupled protein-ligand dynamics in truncated hemoglobin N from atomistic simulations and transition networks

BACKGROUND

The nature of ligand motion in proteins is difficult to characterize directly using experiment. Specifically, it is unclear to what degree these motions are coupled.

METHODS

All-atom simulations are used to sample ligand motion in truncated Hemoglobin N. A transition network analysis including ligand- and protein-degrees of freedom is used to analyze the microscopic dynamics.

RESULTS

Clustering of two different subsets of MD trajectories highlights the importance of a diverse and exhaustive description to define the macrostates for a ligand-migration network. Monte Carlo simulations on the transition matrices from one particular clustering are able to faithfully capture the atomistic simulations. Contrary to clustering by ligand positions only, including a protein degree of freedom yields considerably improved coarse grained dynamics. Analysis with and without imposing detailed balance agree closely which suggests that the underlying atomistic simulations are converged with respect to sampling transitions between neighboring sites.

CONCLUSIONS

Protein and ligand dynamics are not independent from each other and ligand migration through globular proteins is not passive diffusion.

GENERAL SIGNIFICANCE

Transition network analysis is a powerful tool to analyze and characterize the microscopic dynamics in complex systems. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

Hydration free energies of cyanide and hydroxide ions from molecular dynamics simulations with accurate force fields

The evaluation of hydration free energies is a sensitive test to assess force fields used in atomistic simulations. We showed recently that the vibrational relaxation times, 1D- and 2D-infrared spectroscopies for CN(-) in water can be quantitatively described from molecular dynamics (MD) simulations with multipolar force fields and slightly enlarged van der Waals radii for the C- and N-atoms. To validate such an approach, the present work investigates the solvation free energy of cyanide in water using MD simulations with accurate multipolar electrostatics. It is found that larger van der Waals radii are indeed necessary to obtain results close to the experimental values when a multipolar force field is used. For CN(-), the van der Waals ranges refined in our previous work yield hydration free energy between -72.0 and -77.2 kcal mol(-1), which is in excellent agreement with the experimental data. In addition to the cyanide ion, we also study the hydroxide ion to show that the method used here is readily applicable to similar systems. Hydration free energies are found to sensitively depend on the intermolecular interactions, while bonded interactions are less important, as expected. We also investigate in the present work the possibility of applying the multipolar force field in scoring trajectories generated using computationally inexpensive methods, which should be useful in broader parametrization studies with reduced computational resources, as scoring is much faster than the generation of the trajectories.

Diffusion of atomic oxygen relevant to water formation in amorphous interstellar ices

Molecular dynamics (MD) simulations together with accurate physics-based force fields are employed to determine the mobility of atomic oxygen in amorphous ice at low temperatures, characteristic for conditions in interstellar ices. From the simulations it is found that the mobility of atomic oxygen ranges from 60 to 480 Ų ns⁻¹ in amorphous ice at temperatures between 50 and 200 K. Hence, the simulations establish that atomic oxygen is mobile to a certain degree and a chemical mechanism for water formation involving oxygen mobility is a realistic scenario. This is also confirmed by the computed migration barriers for oxygen diffusion by multiple umbrella sampling simulations, which yield barriers for diffusion in the range of 0.7–1.9 kcal mol⁻¹. The physics-based force field – based on a multipolar expansion of the electrostatic interactions – yields more pronounced energetics for oxygen migration pathways compared to the conventional point-charge models employed in typical simulations. Once formed, the computed solvation free energy suggests that atomic oxygen thermodynamically prefers to be localized inside amorphous ice and is available for chemical reaction, which may be relevant to water formation in and on grains.

Vibrational dynamics and solvatochromism of the label SCN in various solvents and hemoglobin by time dependent IR and 2D-IR spectroscopy

The vibrational label SCN is used to report on local structural dynamics in a protein revealing spectral diffusion on a picosecond scale. The SCN spectra are compared to the response of methylthiocyanate in solvents with different polarity and hydrogen-bonding capabilities.