Grand Molecular Dynamics: A Method for Open Systems

We present a new molecular dynamics method for studying the dynamics of open systems. The method couples a classical system to a chemical potential reservior. In the formulation, following the extended system dynamics approach, we introduce a variable, v to represent the coupling to the chemical potential reservoir. The new variable governs the dynamics of the variation of number of particles in the system. The number of particles is determined by taking the integer part of v. The fractional part of the new variable is used to scale the potential energy and the kinetic energy of an additional particle: i.e., we introduce a fractional particle. We give the ansatz Lagrangians and equations of motion for both the isothermal and the adiabatic forms of grand molecular dynamics. The averages calculated over the trajectories generated by these equations of motion represent the classical grand canonical ensemble (μVT) and the constant chemical potential adiabatic ensemble (μVL) averages, respectively. The microcanonical phase space densities of the adiabatic and isothermal forms the molecular dynamics method are shown to be equivalent to adiabatic constant chemical potential ensemble, and grand canonical ensemble partition functions. We also discuss the extension to multi-component systems, molecular fluids, ionic solutions and the problems and solutions associated with the implementation of the method. The statistical expressions for thermodynamic functions such as specific heat; adiabatic bulk modulus, Grüneissen parameter and number fluctuations are derived. These expressions are used to analyse trajectories of constant chemical potential systems.     

Free energy of formation of clusters of sulphuric acid and water molecules determined by guided disassembly

We evaluate the grand potential of a cluster of two molecular species, equivalent to its free energy of formation from a binary vapour phase, using a non-equilibrium molecular dynamics technique where guide particles, each tethered to a molecule by a harmonic force, move apart to disassemble a cluster into its components. The mechanical work performed in an ensemble of trajectories is analysed using the Jarzynski equality to obtain a free energy of disassembly, a contribution to the cluster grand potential. We study clusters of sulphuric acid and water at 300 K, using a classical interaction scheme, and contrast two modes of guided disassembly. In one, the cluster is broken apart through simple pulling by the guide particles, but we find the trajectories tend to be mechanically irreversible. In the second approach, the guide motion and strength of tethering are modified in a way that prises the cluster apart, a procedure that seems more reversible. We construct a surface representing the cluster grand potential, and identify a critical cluster for droplet nucleation under given vapour conditions. We compare the equilibrium populations of clusters with calculations reported by Henschel et al. [J. Phys. Chem. A 2014;118:2599] based on optimised quantum chemical structures.     

Harder,better, faster,stronger: Large-scale QM and QM/MM for predictive modeling in enzymes and proteins

Computational prediction of enzyme mechanism and protein function requires accurate physics-based models and suitable sampling. We discuss recent advances in large-scale quantum mechanical (QM) modeling of biochemical systems that have reduced the cost of high-accuracy models. Tradeoffs between sampling and accuracy have motivated modeling with molecular mechanics (MM) in a multiscale QM/MM or iterative approach. Limitations to both conventional density-functional theory and classical MM force fields remain for describing noncovalent interactions in comparison to experiment or wavefunction theory. Because predictions of enzyme action (i.e. electrostatics), free energy barriers, and mechanisms are sensitive to the protocol and embedding method in QM/MM, convergence tests and systematic methods for quantifying QM-level interactions are a needed, active area of development.     

Fluorescence correlation spectroscopy: past, present, future

In recent years fluorescence correlation spectroscopy (FCS) has become a routine method for determining diffusion coefficients, chemical rate constants, molecular concentrations, fluorescence brightness, triplet state lifetimes, and other molecular parameters. FCS measures the spatial and temporal correlation of individual molecules with themselves and so provides a bridge between classical ensemble and contemporary single-molecule measurements. It also provides information on concentration and molecular number fluctuations for nonlinear reaction systems that complement single-molecule measurements. Typically implemented on a fluorescence microscope, FCS samples femtoliter volumes and so is especially useful for characterizing small dynamic systems such as biological cells. In addition to its practical utility, however, FCS provides a window on mesoscopic systems in which fluctuations from steady states not only provide the basis for the measurement but also can have important consequences for the behavior and evolution of the system. For example, a new and potentially interesting field for FCS studies could be the study of nonequilibrium steady states, especially in living cells.     

Linking Well-Tempered Metadynamics Simulations with Experiments

Linking experiments with the atomistic resolution provided by molecular dynamics simulations can shed light on the structure and dynamics of protein-disordered states. The sampling limitations of classical molecular dynamics can be overcome using metadynamics, which is based on the introduction of a history-dependent bias on a small number of suitably chosen collective variables. Even if such bias distorts the probability distribution of the other degrees of freedom, the equilibrium Boltzmann distribution can be reconstructed using a recently developed reweighting algorithm. Quantitative comparison with experimental data is thus possible. Here we show the potential of this combined approach by characterizing the conformational ensemble explored by a 13-residue helix-forming peptide by means of a well-tempered metadynamics/parallel tempering approach and comparing the reconstructed nuclear magnetic resonance scalar couplings with experimental data.     

Molecular Dynamics with Nonadiabatic Transitions: A Comparison of Methods

Abstract

Surface hopping (SH) and density matrix evolution (DME) methods which simulate the dynamics of quantum systems embedded in a classical environments are compared with exact quantum-dynamical calculations. These methods are applied to study the inelastic collisions of a classical particle with a five-level quantum harmonic oscillator. One-dimensional, two-state models representing electronic transitions are also treated. In addition, the methods are applied to the dynamics of a proton in a bistable potential bilinearly coupled to the bath of classical harmonic oscillators. Vibrational spectra calculated by both methods compare well with each other. The SH results are, in general, closer to the results of a full quantum treatment than the corresponding DME values. The DME method breaks down in the case of extended coupling with reflection at low energies.     

A milestoning study of the kinetics of an allosteric transition: atomically detailed simulations of deoxy Scapharca hemoglobin

Atomically detailed simulations are used to compute the kinetics of the R-to-T transition in deoxy Scapharca hemoglobin. A computational approach called milestoning is utilized that combines 1), an efficient reaction path algorithm; and 2), a "fragment and glue" approach for classical trajectories. Milestoning computes the R-to-T transition kinetics on the microsecond timescale based on atomically detailed trajectories that rarely exceed a nanosecond. Eleven reference hypersurfaces (milestones) are constructed along the reaction coordinate, which is computed with a global path optimization algorithm. Two-hundred classical trajectories are calculated for each of the milestones to collect local distributions of first passage times. These local distributions are used in a non-Markovian theory to compute the overall timescale. Exponential enrichment of reactive trajectories, an important component of the milestoning approach, makes these calculations possible. The overall timescale of the reaction is estimated as 10 +/- 9 micros, in accord with available experimental data. The barrier is not sharp and is spread over four milestones. Even after the most significant structural changes are completed (phenylalanine F4 ring flips), highly collective and activated motions continue. The calculations suggest an additional late free energy barrier.     

Comparison of molecular interactions of Ag2Te and CdTe quantum dots with human serum albumin by spectroscopic approaches

Ag₂Te quantum dots (QDs) have attracted great attention in biological applications due to their superior photoluminescence qualities and good biocompatibility, but their potential biotoxicity at a molecular biology level has been rarely discussed. In order to better understand the basic behavior of Ag₂Te QDs in biological systems and compare their biotoxicity to cadmium‐containing QDs, a series of spectroscopic measurements was applied to reveal the molecular interactions of Ag₂Te QDs and CdTe QDs with human serum albumin (HSA). Ag₂Te QDs and CdTe QDs statically quenched the intrinsic fluorescence of HSA by electrostatic interactions, but Ag₂Te QDs exhibited weaker quenching ability and weaker binding ability compared with CdTe QDs. Electrostatic interactions were the main binding forces and Sudlow's site I was the primary binding site during these binding interactions. Furthermore, micro‐environmental and conformational variations of HSA were induced by their binding interactions with two QDs. Ag₂Te QDs caused less secondary structural and conformational change in HSA, illustrating the lower potential biotoxicity risk of Ag₂Te QDs. Our results systematically indicated the molecular binding mechanism of Ag₂Te QDs with HSA, which provided important information for possible toxicity risk of these cadmium‐free QDs to human health.     

Probing the initial stage of aggregation of the Abeta(10-35)-protein: assessing the propensity for peptide dimerization

Characterization of the early stages of peptide aggregation is of fundamental importance in elucidating the mechanism of the formation of deposits associated with amyloid disease. The initial step in the pathway of aggregation of the Abeta-protein, whose monomeric NMR structure is known, was studied through the simulation of the structure and stability of the peptide dimer in aqueous solution. A protocol based on shape complementarity was used to generate an assortment of possible dimer structures. The structures generated based on shape complementarity were evaluated using rapidly computed estimates of the desolvation and electrostatic interaction energies to identify a putative stable dimer structure. The potential of mean force associated with the dimerization of the peptides in aqueous solution was computed for both the hydrophobic and the electrostatic driven forces using umbrella sampling and classical molecular dynamics simulation at constant temperature and pressure with explicit solvent and periodic boundary conditions. The comparison of the two free energy profiles suggests that the structure of the peptide dimer is determined by the favorable desolvation of the hydrophobic residues at the interface. Molecular dynamics trajectories originating from two putative dimer structures indicate that the peptide dimer is stabilized primarily through hydrophobic interactions, while the conformations of the peptide monomers undergo substantial structural reorganization in the dimerization process. The finding that the phi-dimer may constitute the ensemble of stable Abeta(10-35) dimer has important implications for fibril formation. In particular, the expulsion of water molecules at the interface might be a key event, just as in the oligomerization of Abeta(16-22) fragments. We conjecture that events prior to the nucleation process themselves might involve crossing free energy barriers which depend on the peptide-peptide and peptide-water interactions. Consistent with existing experimental studies, the peptides within the ensemble of aggregated states show no signs of formation of secondary structure.     

Constant electric field simulations of the membrane potential illustrated with simple systems

Advances in modern computational methods and technology make it possible to carry out extensive molecular dynamics simulations of complex membrane proteins based on detailed atomic models. The ultimate goal of such detailed simulations is to produce trajectories in which the behavior of the system is as realistic as possible. A critical aspect that requires consideration in the case of biological membrane systems is the existence of a net electric potential difference across the membrane. For meaningful computations, it is important to have well validated methodologies for incorporating the latter in molecular dynamics simulations. A widely used treatment of the membrane potential in molecular dynamics consists of applying an external uniform electric field E perpendicular to the membrane. The field acts on all charged particles throughout the simulated system, and the resulting applied membrane potential V is equal to the applied electric field times the length of the periodic cell in the direction perpendicular to the membrane. A series of test simulations based on simple membrane-slab models are carried out to clarify the consequences of the applied field. These illustrative tests demonstrate that the constant-field method is a simple and valid approach for accounting for the membrane potential in molecular dynamics studies of biomolecular systems. This article is part of a Special Issue entitled: Membrane protein structure and function.     

Virtual and solution conformations of oligosaccharides

The possibility that observed nuclear Overhauser enhancements and bulk longitudinal relaxation times, parameters measured by 1H NMR and often employed in determining the preferred solution conformation of biologically important molecules, are the result of averaging over many conformational states is quantitatively evaluated. Of particular interest was to ascertain whether certain 1H NMR determined conformations are "virtual" in nature; i.e., the fraction of the population of molecules actually found at any time within the subset of conformational space defined as the "solution conformation" is vanishingly small. A statistical mechanics approach was utilized to calculate an ensemble average relaxation matrix from which (NOE)'s and (T1)'s are calculated. Model glycosidic linkages in four oligosaccharides were studied. The solution conformation at any glycosidic linkage is properly represented by a normalized, Boltzmann distribution of conformers generated from an appropriate potential energy surface. The nature of the resultant population distributions is such that 50% of the molecular population is found within 1% of available microstates, while 99% of the molecular population occupies about 10% of the ensemble microstates, a number roughly equal to that sterically allowed. From this analysis we conclude that in many cases quantitative interpretation of NMR relaxation data, which attempts to define a single set of allowable torsion angle values consistent with the observed data, will lead to solution conformations that are either virtual or reflect torsion angle values possessed by a minority of the molecular population. On the other hand, calculation of ensemble average NMR relaxation data yields values in agreement with experimental results. Observed values of NMR relaxation data are the result of the complex interdependence of the population distribution and NOE (or T1) surfaces in conformational space. In conformational analyses, NMR data can therefore be used to test different population distributions calculated from empirical potential energy functions.     

Heat capacities of supercritical fluids via Grand Canonical ensemble simulations

In this work, suitable mathematical relationships to compute isobaric heat capacities from molecular simulations in the Grand Canonical (GC) ensemble are derived and tested via Monte Carlo methods. Using atomistic classical force fields, the residual isobaric heat capacities of pure carbon dioxide (CO₂) and pure methanol (MeOH) were obtained at supercritical conditions (with critical properties estimated from a finite-size scaling analysis). The total isobaric heat capacity was determined by combining the residual isobaric heat capacity obtained from molecular simulations with the ideal gas contributions obtained from experimental correlations. Isobaric heat capacities generated from both GC and Isothermal–Isobaric ensemble simulations were compared to predictions from accurate equations of state (EOS)s for CO₂ and MeOH at corresponding reduced temperatures and pressures. Isobaric heat capacities calculated from both ensembles were in good agreement with those obtained from the Span and Wagner EOS for CO₂ and the IUPAC EOS for MeOH. For comparable computation times, simulations run in the GC ensemble generate results with significantly lower statistical uncertainty than those run in the Isothermal–Isobaric ensemble.     

Quantum simulation for muoniated and deuterated methyl radicals in implicit water solvent: combined ab initio path integral molecular dynamics and the polarizable continuum model simulation study

We have combined ab initio path integral molecular dynamics (PIMD) simulation and the polarizable continuum model (PCM) method to efficiently incorporate solvent effects into nuclear quantum fluctuation of molecular systems. Our combined ab initio PIMD–PCM simulation was applied to muoniated and deuterated methyl radical immersed in implicit water solvent to gain information on solvent and isotope effects from one simulation run. We found that solvent effects lead to the bond elongation and a decrease in the magnitude of isotropic hyperfine coupling constants. These are consistent with the trends in conventional static calculations and experiments. In addition, the performance of cavity models (universal force field, united atom specified for Kohn–Sham and these hybrid models) and the conservation of the PIMD–PCM Hamiltonian were accessed. We confirmed that solvent effects on nuclear quantum fluctuation are efficiently computed using our combined simulation of quantum solute in implicit solvent.     

Discriminating D1 and D2 agonists with a hydrophobic similarity index

Currently, methods for calculating molecular similarity indices have been developed for comparing steric, charge density, and molecular electrostatic potential (MEP) properties. Much of the existing technology may, however, be applied to the quantitative comparison of molecular hydrophobicities. In this article we present an empirical hydrophobic similarity index. We utilize atomic hydrophobic parameters derived from a quantum mechanical semiempirical wavefunction. Hydrophobicity at points on a grid is computed with a recently introduced “molecular lipophilicity potential”. The overlap of pairs of molecules is calculated with the metric introduced by Carbó. This approach is applied to a case in which steric and electrostatic criteria have already been shown to be inadequate in rationalizing selectivity, namely, requirements for recognition at the dopamine D1 and D2 receptors. We demonstrate that, for a set of dopamine agonists, D1 ligands show higher similarity in this property than D2 analogs. This indicator of similarity is more successful at accounting for D1 selectivity than previous methods.     

Quantum-dynamical picture of a multistep enzymatic process: reaction catalyzed by phospholipase A(2)

A quantum-classical molecular dynamics model (QCMD), applying explicit integration of the time-dependent Schr?dinger equation (QD) and Newtonian equations of motion (MD), is presented. The model is capable of describing quantum dynamical processes in complex biomolecular systems. It has been applied in simulations of a multistep catalytic process carried out by phospholipase A(2) in its active site. The process includes quantum-dynamical proton transfer from a water molecule to histidine localized in the active site, followed by a nucleophilic attack of the resulting OH(-) group on a carbonyl carbon atom of a phospholipid substrate, leading to cleavage of an adjacent ester bond. The process has been simulated using a parallel version of the QCMD code. The potential energy function for the active site is computed using an approximate valence bond (AVB) method. The dynamics of the key proton is described either by QD or classical MD. The coupling between the quantum proton and the classical atoms is accomplished via Hellmann-Feynman forces, as well as the time dependence of the potential energy function in the Schr?dinger equation (QCMD/AVB model). Analysis of the simulation results with an Advanced Visualization System revealed a correlated rather than a stepwise picture of the enzymatic process. It is shown that an sp(2)--> sp(3) configurational change at the substrate carbonyl carbon is mostly responsible for triggering the activation process.     

Receptors and functional linkage in membrane permeability: A quantum mechanical model

During functional linkage, ligand receptors are coupled to other receptors and to the cell's metabolic-transport apparatus. The linkage guides the cellular processing of matter, energy and information. Previous conceptions of functional linkage have used the ideas of classical physics appropriate to macroscopic objects. This study presents an initial quantum mechanical model of functional linkage in the case of ligands moving through lipid bilayers and hydrophilic transmembrane channels (‘pores’) of molecular dimensions. On the basis of permeability data, energy surfaces consisting of piecewise-constant potential regions are used to model the lipid bilayers and transmembrane channels. The centre-of-mass wavefunction for a ligand on such energy surfaces is analysed and the permeability coefficients calculated from the wavefunction's transmission characteristics. It is found that quasi-bound states in the several ligand-binding regions of a bilayer or pore system can functionally link to facilitate the passage of the molecule across the permeability barrier. Appearance of the linkage is a sensitive function of the ligand's energy. If the centre-of-mass energies are distributed as in a thermalized fluid, the flux via the quantum functional linkage can equal or exceed that of a classical flux for proton transport through rigid pores in which the intrasite barriers are relatively high (0.25–1 eV) and narrow (0.1–1 Å). The functional linkage plays a less important role in bilayer (rather than pore) energy surfaces and at higher molecular weights. If the ligand-receptor interaction is accompanied by energy transfer to or from ligands, the flux via the quantum functional linkage can equal or exceed the classically expected flux at all relevant ligand molecular weights. These findings are discussed in relation to earlier work and the limitations of the model emphasized.     

Abstracts: Albany 2007, The 15th Conversation

Abstract

Forty nine molecular dynamics simulations of unfolding trajectories of the segment B1 of streptococcal protein G (GB1) provide a direct demonstration of the diversity of unfolding pathway and give a statistically utmost unfolding pathway under the physical property space. Twelve physical properties of the protein were chosen to construct a 12-dimensional property space. Then the 12-dimentional property space was reduced to a 3-dimentional principle component property space. Under the property space, the multiple unfolding trajectories look like “trees”, which have some common characters. The “root of the tree” corresponds to the native state, the “bole” homologizes the partially unfolded conformations, and the “crown” is in correspondence to the unfolded state. These unfolding trajectories can be divided into three types. The first one has the characters of straight “bole” and “crown” corresponding to a fast two-state unfolding pathway of GB1. The second one has the character of “the standstill in the middle tree bole”, which may correspond to a three-state unfolding pathway. The third one has the character of “the circuitous bole” corresponding to a slow two-state unfolding pathway. The fast two-state unfolding pathway is a statistically utmost unfolding pathway or preferred pathway of GB1, which occupies 53% of 49 unfolding trajectories. In the property space all the unfolding trajectories construct a thermal unfolding pathway ensemble of GB1. The unfolding pathway ensemble resembles a funnel that is gradually emanative from the native state ensemble to the unfolded state ensemble. In the property space, the thermal unfolded state distribution looks like electronic cloud in quantum mechanics. The unfolded states of the independent unfolding simulation trajectories have substantial overlaps, indicating that the thermal unfolded states are confined by the physical property values, and the number of protein unfolded state are much less than that was believed before.     

Analytical expressions on harmonic oscillators for variational path integrals

In this study, we derive analytical expressions for one-dimensional harmonic oscillators for variational path integrals (VPIs). A Gaussian-type trial wavefunction is adopted. Total and potential energies of the system are analytically expressed both for the continuous time and an approximate discretised VPIs. Obtained expressions are numerically verified using molecular dynamics calculations. Convergence properties regarding the projection time and the Trotter number are discussed.