Advanced Potential Energy Surfaces for Molecular Simulation

Mechanism of thermal oxidation into volatile compounds from (E)-4-decenal: A density functional theory study

Unsaturated aliphatic aldehyde oxidation plays a significant role in the deep oxidation of fatty acids to produce volatile chemicals. Exposing the oxidation process of unsaturated aliphatic aldehydes is crucial to completely comprehend how food flavor forms. In this study, thermal desorption cryo-trapping in conjunction with gas chromatography-mass spectrometry was used to examine the volatile profile of (E)-4-decenal during heating, and 32 volatile compounds in all were detected and identified. Meanwhile, density functional theory (DFT) calculations were used, and 43 reactions were obtained in the 24 pathways, which were summarized into the peroxide reaction mechanism (ROOH), the peroxyl radical reaction mechanism (ROO·) and the alkoxy radical reaction mechanism (RO·). Moreover, the priority of these three oxidative mechanisms was the RO· mechanism > ROOH mechanism > ROO· mechanism. Furthermore, the DFT results and experimental results agreed well, and the oxidative mechanism of (E)-4-decenal was finally illuminated.

The Power of Molecular Dynamics Simulations and Their Applications to Discover Cysteine Protease Inhibitors

A large family of enzymes with the function of hydrolyzing peptide bonds, called peptidases or cysteine proteases (CPs), are divided into three categories according to the peptide chain involved. CPs catalyze the hydrolysis of amide, ester, thiol ester, and thioester peptide bonds. They can be divided into several groups, such as papain-like (CA), viral chymotrypsin-like CPs (CB), papainlike endopeptidases of RNA viruses (CC), legumain-type caspases (CD), and showing active residues of His, Glu/Asp, Gln, Cys (CE). The catalytic mechanism of CPs is the essential cysteine residue present in the active site. These mechanisms are often studied through computational methods that provide new information about the catalytic mechanism and identify inhibitors. The role of computational methods during drug design and development stages is increasing. Methods in Computer-Aided Drug Design (CADD) accelerate the discovery process, increase the chances of selecting more promising molecules for experimental studies, and can identify critical mechanisms involved in the pathophysiology and molecular pathways of action. Molecular dynamics (MD) simulations are essential in any drug discovery program due to their high capacity for simulating a physiological environment capable of unveiling significant inhibition mechanisms of new compounds against target proteins, especially CPs. Here, a brief approach will be shown on MD simulations and how the studies were applied to identify inhibitors or critical information against cysteine protease from several microorganisms, such as Trypanosoma cruzi (cruzain), Trypanosoma brucei (rhodesain), Plasmodium spp. (falcipain), and SARS-CoV-2 (Mpro). We hope the readers will gain new insights and use our study as a guide for potential compound identifications using MD simulations.

Data scheme and data format for transferable force fields for molecular simulation

A generalized data scheme for transferable classical force fields used in molecular simulations, i.e. molecular dynamics and Monte Carlo simulation, is presented. The data scheme is implemented in an SQL-based data format. The data scheme and data format is machine readable, re-usable, and interoperable. A transferable force field is a chemical construction plan specifying intermolecular and intramolecular interactions between different types of atoms or different chemical groups and can be used for building a model for a given component. The data scheme proposed in this work (named TUK-FFDat) formalizes digitally these chemical construction plans, i.e. transferable force fields. It can be applied to all-atom as well as united-atom transferable force fields. The general applicability of the data scheme is demonstrated for different types of force fields (TraPPE, OPLS-AA, and Potoff). Furthermore, conversion tools for translating the data scheme between .xls spread sheet format and the SQL-based data format are provided. The data format can readily be integrated in existing workflows, simulation engines, and force field databases as well as for linking such.

Concatenation of molecular docking and dynamics simulation of human papillomavirus type 16 E7 oncoprotein targeted ligands: In quest of cervical cancer's treatment

The Human papillomaviruses type 16 E7 oncoprotein is a 98-amino-acid, 11-kilodalton acidic oncoprotein with three conserved portions. Due to its interaction with the pRb-E2F complex, CKII, CKI (mostly p21), and even HDAC1, it possesses strong transformative and carcinogenic qualities that inhibit normal differentiation and cell cycle regulation. Here, we target the E7 oncoprotein using two prior research active compounds: asarinin and thiazolo[3,2-a]benzimidazole-3(2H)-one,2-(2-fluorobenzylideno)-7,8-dimethyl (thiazolo), and valproic acid as a control. We are performing molecular docking followed by molecular dynamic analysis. By acting as competitive inhibitors in the binding site, it was hypothesized that both drugs would inhibit E7-mediated pRb degradation and E7-mediated p21 degradation, resulting in decreased cell cycle progression, immortalization, and proliferation. In addition, we expect that the direct inhibitory action of valproic acid in E7 will target the CKII-mediated phosphorylation pathway necessary for destabilizing p130 and pRb. According to the results of the dynamic simulation, stable interactions exist between every compound. Despite the instability of E7 protein, stability results indicate that both natural chemicals are preferable, with thiazolo outperforming valproic acid.

Force Decomposition Analysis: A Method to Decompose Intermolecular Forces into Physically Relevant Component Contributions

Computational quantum chemistry can be more than just numerical experiments when methods are specifically adapted to investigate chemical concepts. One important example is the development of energy decomposition analysis (EDA) to reveal the physical driving forces behind intermolecular interactions. In EDA, typically the interaction energy from a good-quality density functional theory (DFT) calculation is decomposed into multiple additive components that unveil permanent and induced electrostatics, Pauli repulsion, dispersion, and charge-transfer contributions to noncovalent interactions. Herein, we formulate, implement, and investigate decomposing the forces associated with intermolecular interactions into the same components. The resulting force decomposition analysis (FDA) is potentially useful as a complement to the EDA to understand chemistry, while also providing far more information than an EDA for data analysis purposes such as training physics-based force fields. We apply the FDA based on absolutely localized molecular orbitals (ALMOs) to analyze interactions of water with sodium and chloride ions as well as in the water dimer. We also analyze the forces responsible for geometric changes in carbon dioxide upon adsorption onto (and activation by) gold and silver anions. We also investigate how the force components of an EDA-based force field for water clusters, namely MB-UCB, compare to those from force decomposition analysis.

Learning pair potentials using differentiable simulations

Learning pair interactions from experimental or simulation data is of great interest for molecular simulations. We propose a general stochastic method for learning pair interactions from data using differentiable simulations (DiffSim). DiffSim defines a loss function based on structural observables, such as the radial distribution function, through molecular dynamics (MD) simulations. The interaction potentials are then learned directly by stochastic gradient descent, using backpropagation to calculate the gradient of the structural loss metric with respect to the interaction potential through the MD simulation. This gradient-based method is flexible and can be configured to simulate and optimize multiple systems simultaneously. For example, it is possible to simultaneously learn potentials for different temperatures or for different compositions. We demonstrate the approach by recovering simple pair potentials, such as Lennard-Jones systems, from radial distribution functions. We find that DiffSim can be used to probe a wider functional space of pair potentials compared with traditional methods like iterative Boltzmann inversion. We show that our methods can be used to simultaneously fit potentials for simulations at different compositions and temperatures to improve the transferability of the learned potentials.

Non-covalent interactions from a Quantum Chemical Topology perspective

About half a century after its little-known beginnings, the quantum topological approach called QTAIM has grown into a widespread, but still not mainstream, methodology of interpretational quantum chemistry. Although often confused in textbooks with yet another population analysis, be it perhaps an elegant but somewhat esoteric one, QTAIM has been enriched with about a dozen other research areas sharing its main mathematical language, such as Interacting Quantum Atoms (IQA) or Electron Localisation Function (ELF), to form an overarching approach called Quantum Chemical Topology (QCT). Instead of reviewing the latter's role in understanding non-covalent interactions, we propose a number of ideas emerging from the full consequences of the space-filling nature of topological atoms, and discuss how they (will) impact on interatomic interactions, including non-covalent ones. The architecture of a force field called FFLUX, which is based on these ideas, is outlined. A new method called Relative Energy Gradient (REG) is put forward, which is able, by computation, to detect which fragments of a given molecular assembly govern the energetic behaviour of this whole assembly. This method can offer insight into the typical balance of competing atomic energies both in covalent and non-covalent case studies. A brief discussion on so-called bond critical points is given, highlighting concerns about their meaning, mainly in the arena of non-covalent interactions.

Two Liberibacter Proteins Combine to Suppress Critical Innate Immune Defenses in Citrus

We adopted a systems-based approach to determine the role of two Candidatus Liberibacter asiaticus (CLas) proteins, LasP ₂₃₅ and Effector 3, in Huanglongbing (HLB) pathogenesis. While a published work suggests the involvement of these CLas proteins HLB pathogenesis, the exact structure-based mechanism of their action has not been elucidated. We conducted the following experiments to determine the structure-based mechanisms of action. First, we immunoprecipitated the interacting citrus protein partners of LasP ₂₃₅ and Effector 3 from the healthy and CLas-infected Hamlin extracts and identified them by Liquid Chromatography with tandem mass spectrometry (LC-MS/MS). Second, we performed a split green fluorescent protein (GFP) assay in tobacco to validate that the interactions observed in vitro are also retained in planta. The notable in planta citrus targets of LasP ₂₃₅ and Effector 3 include citrus innate immune proteins. Third, in vitro and in planta studies were performed to show that LasP ₂₃₅ and Effector 3 interact with and inhibit the functions of multiple citrus proteins belonging to the innate immune pathways. These inhibitory interactions led to a high level of reactive oxygen species, blocking of bactericidal lipid transfer protein (LTP), and induction of premature programed cell death (PCD), all of which are beneficial to CLas lifecycle and HLB pathogenesis. Finally, we performed molecular dynamics simulations to visualize the interactions of LasP ₂₃₅ and Effector 3, respectively, with LTP and Kunitz protease inhibitor. This led to the design of an LTP mimic, which sequestered and blocked LasP ₂₃₅ and rescued the bactericidal activity of LTP thereby proving that LasP ₂₃₅ , indeed, participates in HLB pathogenesis.

Progress in understanding crystallisation: a personal perspective

After this discussion meeting, most participants felt that we do not understand crystallisation. However, in the 1980s I believe that most scientists would have considered that crystallisation was adequately understood....

Protein-Ligand Binding Molecular Details Revealed by Terahertz Optical Kerr Spectroscopy: A Simulation Study

Picosecond fast motions and their involvement in the biochemical processes such as protein-ligand binding has engaged significant attention. Terahertz optical Kerr spectroscopy (OKE) has the superior potential to probe these fast motions directly. Application of OKE in protein-ligand binding study is, however, limited by the difficulty of quantitative atomistic interpretation, and the calculation of Kerr spectrum for entire solvated protein complex was considered not yet feasible, due to the lack of one consistent polarizable model for both configuration sampling and polarizability calculation. Here, we analyzed the biochemical relevance of OKE to the lysozyme-triacetylchitotriose binding based on the first OKE simulation using one consistent Drude polarizable model. An analytical multipole and induced dipole scheme was employed to calculate the off-diagonal Drude polarizability more efficiently and accurately. Further theoretical analysis revealed how the subtle twisting and stiffening of aromatic protein residues' spatial arrangement as well as the confinement of small water clusters between ligand and protein cavity due to the ligand binding can be examined using Kerr spectroscopy. Comparison between the signals of bound complex and that of uncorrelated protein/ligand demonstrated that binding action alone has reflection in the OKE spectrum. Our study indicated OKE as a powerful terahertz probe for protein-ligand binding chemistry and dynamics.

Reduced-dimensional surface hopping with offline-online computations

Molecular dynamics simulations often classically evolve the nuclear geometry on adiabatic potential energy surfaces (PESs), punctuated by random hops between energy levels in regions of strong coupling, in an algorithm known as surface hopping. However, the computational expense of integrating the geometry on a full-dimensional PES and computing the required couplings can quickly become prohibitive as the number of atoms increases. In this work, we describe a method for surface hopping that uses only important reaction coordinates, performs all expensive evaluations of the true PESs and couplings only once before simulating dynamics (offline), and then queries the stored values during the surface hopping simulation (online). Our Python codes are freely available on GitHub. Using photodissociation of azomethane as a test case, this method is able to reproduce experimental results that have thus far eluded ab initio surface hopping studies.

Efficient Approximation of Potential Energy Surfaces with Mixed-Basis Interpolation

The potential energy surface (PES) describes the energy of a chemical system as a function of its geometry and is a fundamental concept in modern chemistry. A PES provides much useful information about the system, including the structures and energies of various stationary points, such as stable conformers (local minima) and transition states (first-order saddle points) connected by a minimum-energy path. Our group has previously produced surrogate reduced-dimensional PESs using sparse interpolation along chemically significant reaction coordinates, such as bond lengths, bond angles, and torsion angles. These surrogates used a single interpolation basis, either polynomials or trigonometric functions, in every dimension. However, relevant molecular dynamics (MD) simulations often involve some combination of both periodic and nonperiodic coordinates. Using a trigonometric basis on nonperiodic coordinates, such as bond lengths, leads to inaccuracies near the domain boundary. Conversely, polynomial interpolation on the periodic coordinates does not enforce the periodicity of the surrogate PES gradient, leading to nonconservation of total energy even in a microcanonical ensemble. In this work, we present an interpolation method that uses trigonometric interpolation on the periodic reaction coordinates and polynomial interpolation on the nonperiodic coordinates. We apply this method to MD simulations of possible isomerization pathways of azomethane between cis and trans conformers. This method is the only known interpolative method that appropriately conserves total energy in systems with both periodic and nonperiodic reaction coordinates. In addition, compared to all-polynomial interpolation, the mixed basis requires fewer electronic structure calculations to obtain a given level of accuracy, is an order of magnitude faster, and is freely available on GitHub.

Effective Two-Body Interactions

Cooperative or nonadditive effects contribute to the pairwise noncovalent interaction of two molecules in a cluster or the condensed phase in ways that depend on the specific arrangements and interactions of the other surrounding molecules that constitute their environment. General expressions for an effective two-body interaction are presented, which are correct to increasing orders in the many-body expansion. The simplest result, correct through third order, requires only seven individual calculations, in contrast to a linear number of three-body contributions. Two applications are presented. First, an error analysis is performed on a model (H₂O)₈ cluster which completes the first solvation shell of a central water-water hydrogen bond. Energy decomposition analysis is performed to show that the largest effects of cooperativity on the central hydrogen bond arise from electrical polarization. Second, the nature of cooperative effects on proton transfer in an HCl + (H₂O)₄ cluster is characterized.

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.

Extension of an Atom-Atom Dispersion Function to Halogen Bonds and Its Use for Rational Design of Drugs and Biocatalysts

A dispersion function Das in the form of a damped atom-atom asymptotic expansion fitted to ab initio dispersion energies from symmetry-adapted perturbation theory was improved and extended to systems containing heavier halogen atoms. To illustrate its performance, the revised Das function was implemented in the multipole first-order electrostatic and second-order dispersion (MED) scoring model. The extension has allowed applications to a much larger set of biocomplexes than it was possible with the original Das. A reasonable correlation between MED and experimentally determined inhibitory activities was achieved in a number of test cases, including structures featuring nonphysically shortened intermonomer distances, which constitute a particular challenge for binding strength predictions. Since the MED model is also computationally efficient, it can be used for reliable and rapid assessment of the ligand affinity or multidimensional scanning of amino acid side-chain conformations in the process of rational design of novel drugs or biocatalysts.

From Intermolecular Interaction Energies and Observable Shifts to Component Contributions and Back Again: A Tale of Variational Energy Decomposition Analysis

Quantum chemistry in the form of density functional theory (DFT) calculations is a powerful numerical experiment for predicting intermolecular interaction energies. However, no chemical insight is gained in this way beyond predictions of observables. Energy decomposition analysis (EDA) can quantitatively bridge this gap by providing values for the chemical drivers of the interactions, such as permanent electrostatics, Pauli repulsion, dispersion, and charge transfer. These energetic contributions are identified by performing DFT calculations with constraints that disable components of the interaction. This review describes the second-generation version of the absolutely localized molecular orbital EDA (ALMO-EDA-II). The effects of different physical contributions on changes in observables such as structure and vibrational frequencies upon complex formation are characterized via the adiabatic EDA. Example applications include red- versus blue-shifting hydrogen bonds; the bonding and frequency shifts of CO, N₂, and BF bound to a [Ru(II)(NH₃)₅]^{2 +} moiety; and the nature of the strongly bound complexes between pyridine and the benzene and naphthalene radical cations. Additionally, the use of ALMO-EDA-II to benchmark and guide the development of advanced force fields for molecular simulation is illustrated with the recent, very promising, MB-UCB potential.

On the faithfulness of molecular mechanics representations of proteins towards quantum-mechanical energy surfaces

Force fields based on molecular mechanics (MM) are the main computational tool to study the relationship between protein structure and function at the molecular level. To validate the quality of such force fields, high-level quantum-mechanical (QM) data are employed to test their capability to reproduce the features of all major conformational substates of a series of blocked amino acids. The phase-space overlap between MM and QM is quantified in terms of the average structural reorganization energies over all energy minima. Here, the structural reorganization energy is the MM potential-energy difference between the structure of the respective QM energy minimum and the structure of the closest MM energy minimum. Thus, it serves as a measure for the relative probability of visiting the QM minimum during an MM simulation. We evaluate variants of the AMBER, CHARMM, GROMOS and OPLS biomolecular force fields. In addition, the two blocked amino acids alanine and serine are used to demonstrate the dependence of the measured agreement on the QM method, the phase, and the conformational preferences. Blocked serine serves as an example to discuss possible improvements of the force fields, such as including polarization with Drude particles, or using tailored force fields. The results show that none of the evaluated force fields satisfactorily reproduces all energy minima. By decomposing the average structural reorganization energies in terms of individual energy terms, we can further assess the individual weaknesses of the parametrization strategies of each force field. The dominant problem for most force fields appears to be the van der Waals parameters, followed to a lesser degree by dihedral and bonded terms. Our results show that performing a simple QM energy optimization from an MM-optimized structure can be a first test of the validity of a force field for a particular target molecule.

Analysis of Photosynthetic Systems and Their Applications with Mathematical and Computational Models

In biological and life science applications, photosynthesis is an important process that involves the absorption and transformation of sunlight into chemical energy. During the photosynthesis process, the light photons are captured by the green chlorophyll pigments in their photosynthetic antennae and further funneled to the reaction center. One of the most important light harvesting complexes that are highly important in the study of photosynthesis is the membrane-attached Fenna–Matthews–Olson (FMO) complex found in the green sulfur bacteria. In this review, we discuss the mathematical formulations and computational modeling of some of the light harvesting complexes including FMO. The most recent research developments in the photosynthetic light harvesting complexes are thoroughly discussed. The theoretical background related to the spectral density, quantum coherence and density functional theory has been elaborated. Furthermore, details about the transfer and excitation of energy in different sites of the FMO complex along with other vital photosynthetic light harvesting complexes have also been provided. Finally, we conclude this review by providing the current and potential applications in environmental science, energy, health and medicine, where such mathematical and computational studies of the photosynthesis and the light harvesting complexes can be readily integrated.

Secondary structure analysis of peptides with relevance to iron-sulfur cluster nesting

Peptides coordinated to iron-sulfur clusters, referred to as maquettes, represent a synthetic strategy for constructing biomimetic models of iron-sulfur metalloproteins. These maquettes have been successfully employed as building blocks of engineered heme-containing proteins with electron-transfer functionality; however, they have yet to be explored in reactivity studies. The concept of iron-sulfur nesting in peptides is a leading hypothesis in Origins-of-Life research as a plausible path to bridge the discontinuity between prebiotic chemical transformations and extant enzyme catalysis. Based on past biomimetic and biochemical research, we put forward a mechanism of maquette reconstitution that guides our development of computational tools and methodologies. In this study, we examined a key feature of the first stage of maquette formation, which is the secondary structure of aqueous peptide models using molecular dynamics simulations based on the AMBER99SB empirical force field. We compared and contrasted S…S distances, [2Fe-2S] and [4Fe-4S] nests, and peptide conformations via Ramachandran plots for dissolved Cys and Gly amino acids, the CGGCGGC 7-mer, and the GGCGGGCGGCGGW 16-mer peptide. Analytical tools were developed for following the evolution of secondary structural features related to [Fe-S] cluster nesting along 100 ns trajectories. Simulations demonstrated the omnipresence of peptide nests for preformed [2Fe-2S] clusters; however, [4Fe-4S] cluster nests were observed only for the 16-mer peptide with lifetimes of a few nanoseconds. The origin of the [4Fe-4S] nest and its stability was linked to a "kinked-ribbon" peptide conformation. Our computational approach lays the foundation for transitioning into subsequent stages of maquette reconstitution, those being the formation of iron ion/iron-sulfur coordinated peptides. © 2018 Wiley Periodicals, Inc.

A FFLUX Water Model: Flexible, Polarizable and with a Multipolar Description of Electrostatics

Key to progress in molecular simulation is the development of advanced models that go beyond the limitations of traditional force fields that employ a fixed, point charge-based description of electrostatics. Taking water as an example system, the FFLUX framework is shown capable of producing models that are flexible, polarizable and have a multipolar description of the electrostatics. The kriging machine-learning methods used in FFLUX are able to reproduce the intramolecular potential energy surface and multipole moments of a single water molecule with chemical accuracy using as few as 50 training configurations. Molecular dynamics simulations of water clusters (25-216 molecules) using the new FFLUX model reveal that incorporating charge-quadrupole, dipole-dipole, and quadrupole-charge interactions into the description of the electrostatics results in significant changes to the intermolecular structuring of the water molecules. © 2019 The Authors. Journal of Computational Chemistry published by Wiley Periodicals, Inc.

Periodic Coulomb Tree Method: An Alternative to Parallel Particle Mesh Ewald

Particle mesh Ewald (PME) is generally the method of choice for handling electrostatics in simulations with periodic boundary conditions. The excellent efficiency of PME on low processor counts is largely due to the use of the fast Fourier transform (FFT). However, due to the FFT's high communication cost, PME scales poorly in parallel. We develop a periodic Coulomb tree (PCT) method for electrostatic interactions in periodic boundary conditions as an alternative to PME in parallel simulations. We verify the accuracy of PCT by comparison of structural and dynamical properties of three different systems obtained via MD simulations using PME and PCT and provide parallel timing comparisons of the two methods on up to 1024 cores.

Approximating Periodic Potential Energy Surfaces with Sparse Trigonometric Interpolation

The potential energy surface (PES) describes the energy of a chemical system as a function of its geometry and is a fundamental concept in computational chemistry. A PES provides much useful information about the system, including the structures and energies of various stationary points, such as local minima, maxima, and transition states. Construction of full-dimensional PESs for molecules with more than 10 atoms is computationally expensive and often not feasible. Previous work in our group used sparse interpolation with polynomial basis functions to construct a surrogate reduced-dimensional PESs along chemically significant reaction coordinates, such as bond lengths, bond angles, and torsion angles. However, polynomial interpolation does not preserve the periodicity of the PES gradient with respect to angular components of geometry, such as torsion angles, which can lead to nonphysical phenomena. In this work, we construct a surrogate PES using trigonometric basis functions, for a system where the selected reaction coordinates all correspond to the torsion angles, resulting in a periodically repeating PES. We find that a trigonometric interpolation basis not only guarantees periodicity of the gradient but also results in slightly lower approximation error than polynomial interpolation.

Natural selection based on coordination chemistry: computational assessment of [4Fe-4S]-maquettes with non-coded amino acids

Cysteine is the only coded amino acid in biology that contains a thiol functional group. Deprotonated thiolate is essential for anchoring iron-sulfur ([Fe-S]) clusters, as prosthetic groups to the protein matrix. [Fe-S] metalloproteins and metalloenzymes are involved in biological electron transfer, radical chemistry, small molecule activation and signalling. These are key metabolic and regulatory processes that would likely have been present in the earliest organisms. In the context of emergence of life theories, the selection and evolution of the cysteine-specific R-CH₂-SH side chain is a fascinating question to confront. We undertook a computational [4Fe-4S]-maquette modelling approach to evaluate how side chain length can influence [Fe-S] cluster binding and stability in short 7-mer and long 16-mer peptides, which contained either thioglycine, cysteine or homocysteine. Force field-based molecular dynamics simulations for [4Fe-4S] cluster nest formation were supplemented with density functional theory calculations of a ligand-exchange reaction between a preassembled cluster and the peptide. Secondary structure analysis revealed that peptides with cysteine are found with greater frequency nested to bind preformed [4Fe-4S] clusters. Additionally, the presence of the single methylene group in cysteine ligands mitigates the steric bulk, maintains the H-bonding and dipole network, and provides covalent Fe-S(thiolate) bonds that together create the optimal electronic and geometric structural conditions for [4Fe-4S] cluster binding compared to thioglycine or homocysteine ligands. Our theoretical work forms an experimentally testable hypothesis of the natural selection of cysteine through coordination chemistry.

Molecular Dynamics Simulations on Relaxed Reduced-Dimensional Potential Energy Surfaces

Molecular dynamics (MD) simulations with full-dimensional potential energy surfaces (PESs) obtained from high-level ab initio calculations are frequently used to model reaction dynamics of small molecules (i.e., molecules with up to 10 atoms). Construction of full-dimensional PESs for larger molecules is, however, not feasible since the number of ab initio calculations required grows rapidly with the increase of dimension. Only a small number of coordinates are often essential for describing the reactivity of even very large systems, and reduced-dimensional PESs with these coordinates can be built for reaction dynamics studies. While analytical methods based on transition-state theory framework are well established for analyzing the reduced-dimensional PESs, MD simulation algorithms capable of generating trajectories on such surfaces are more rare. In this work, we present a new MD implementation that utilizes the relaxed reduced-dimensional PES for standard microcanonical (NVE) and canonical (NVT) MD simulations. The method is applied to the pyramidal inversion of a NH₃ molecule. The results from the MD simulations on a reduced, three-dimensional PES are validated against the ab initio MD simulations, as well as MD simulations on full-dimensional PES and experimental data.

Combining Iteration-Free Polarization with Large Time Step Stochastic-Isokinetic Integration

In order to accelerate molecular dynamics simulations using polarizable force fields, we combine a new extended Lagrangian approach that eliminates the self-consistent field step (iEL/0-SCF) with a stochastic integration scheme that allows for a long time step using a multiple time scale algorithm (SIN(R)). We consider different algorithms for the combined scheme that places different components of the nonbonded forces into different time scales, as well as splitting individual nonbonded forces across time scales, to demonstrate that the combined method works well for bulk water as well as for a concentrated salt solution, aqueous peptide, and solvated protein. Depending on system and desired accuracy, the iEL/0-SCF and SIN(R) combination yields lower bound computational speed-ups of ∼6-8 relative to a molecular dynamics Verlet integration using a standard SCF solver implemented in the reference program TINKER 8.1. The combined approach embodies a significant advance for equilibrium simulations in the canonical ensemble of many-body potential energy surfaces for condensed phase systems with speed-ups that exceed what is possible by either method alone.

Deriving Force-Field Parameters from First Principles Using a Polarizable and Higher Order Dispersion Model

In this work we propose a strategy based on quantum mechanical (QM) calculations to parametrize a polarizable force field for use in molecular dynamics (MD) simulations. We investigate the use of multiple atoms-in-molecules (AIM) strategies to partition QM determined molecular electron densities into atomic subregions. The partitioned atomic densities are subsequently used to compute atomic dispersion coefficients from effective exchange-hole-dipole moment (XDM) calculations. In order to derive values for the repulsive van der Waals parameters from first principles, we use a simple volume relation to scale effective atomic radii. Explicit inclusion of higher order dispersion coefficients was tested for a series of alkanes, and we show that combining C6 and C8 attractive terms together with a C11 repulsive potential yields satisfying models when used in combination with our van der Waals parameters and electrostatic and bonded parameters as directly obtained from quantum calculations as well. This result highlights that explicit inclusion of higher order dispersion terms could be viable in simulation, and it suggests that currently available QM analysis methods allow for first-principles parametrization of molecular mechanics models.

Study of the Wilcox torsion balance in solution for a Tröger's base derivative with hexyl-and heptyl substituents using a combined molecular mechanics and quantum chemistry approach

The folding equilibrium of the Wilcox torsion balance in solution has been studied using a molecular mechanics method for sampling the conformational space and semi-empirical and density-functional quantum chemistry methods for characterizing the relative stabilities of various solute-solvent clusters extracted with the aid of the MD-quench technique from the different simulations that were performed. The role of the solvent environment has been analyzed by choosing four solvents of different polarities, namely water, acetone, tetrachloromethane, and n-hexane. In all cases, it is found that the attractive intramolecular interactions in folded conformations are strongly compensated by the increase of the solute-solvent interaction energies when the molecule unfolds. The latter can be well explained by the larger number of solvent molecules that can bind to the Wilcox molecule when in an unfolded conformation. The results of this work therefore support the experimental results of Yang et al. (Nature Chem 5:1006, 2013) that the folding free energy of the Wilcox balance is strongly reduced in solution as compared to the gas phase.

Toward Learned Chemical Perception of Force Field Typing Rules

Molecular mechanics force fields define how the energy and forces in a molecular system are computed from its atomic positions, thus enabling the study of such systems through computational methods like molecular dynamics and Monte Carlo simulations. Despite progress toward automated force field parametrization, considerable human expertise is required to develop or extend force fields. In particular, human input has long been required to define atom types, which encode chemically unique environments that determine which parameters will be assigned. However, relying on humans to establish atom types is suboptimal. Human-created atom types are often developed without statistical justification, leading to over- or under-fitting of data. Human-created types are also difficult to extend in a systematic and consistent manner when new chemistries must be modeled or new data becomes available. Finally, human effort is not scalable when force fields must be generated for new (bio)polymers, compound classes, or materials. To remedy these deficiencies, our long-term goal is to replace human specification of atom types with an automated approach, based on rigorous statistics and driven by experimental and/or quantum chemical reference data. In this work, we describe novel methods that automate the discovery of appropriate chemical perception: SMARTY allows for the creation of atom types, while SMIRKY goes further by automating the creation of fragment (nonbonded, bonds, angles, and torsions) types. These approaches enable the creation of move sets in atom or fragment type space, which are used within a Monte Carlo optimization approach. We demonstrate the power of these new methods by automating the rediscovery of human defined atom types (SMARTY) or fragment types (SMIRKY) in existing small molecule force fields. We assess these approaches using several molecular data sets, including one which covers a diverse subset of the DrugBank database.

Inertial extended-Lagrangian scheme for solving charge equilibration models

The inertial EL/SCF method is developed to solve charge equilibration models for molecular dynamics, reducing the number of SCFs by 50–80% at each time step when compared to a conjugate gradient SCF solver and tested on diverse reactive systems.

Investigating coordination flexibility of glycerophosphodiesterase (GpdQ) through interactions with mono-, di-, and triphosphoester (NPP, BNPP, GPE, and paraoxon) substrates

Interactions of the catalytically active binuclear form of glycerophosphodiesterase (GpdQ) with chemically diverse substrates, i.e. phosphomono-, phosphodi-, and phosphotriester have been investigated using molecular dynamics (MD) simulations.

β-Lactoglobulin as a potential carrier for bioactive molecules

In this study, the interaction and binding behavior of anesthetic tetracaine (TET) with bovine β-lactoglobulin (LGB) isoform A and a mixture of isoforms A and B were investigated under varying environmental conditions (pH, ionic strength, concentration, LGB-TET complex molar ratio). A wide range of analytical techniques (dynamic light scattering (DLS), electrophoretic mobility, UV-Vis spectroscopy, circular dichroism (CD), quartz crystal microbalance (QCM-D) were used to analyze the physicochemical properties of the complexes in bulk solution and on the surface of gold. The experiments revealed that TET, which is amphiphilic, could bind with LGB not only in the β-barrel but also onto the surface. The zeta potential of the LGB becomes more positively charged upon interaction with TET due to electrostatic interaction of the amino group present in the TET structure. Changes in the zeta potential values are mostly visible above pH 6 for all tested systems. CD spectra show that interaction with the ligand does not change the secondary structure of the protein. The physicochemical properties of LGB-TET complex were examined in an adsorbed state on a gold surface using the QCM-D method. Additionally, molecular docking was used to evaluate the most likely binding site for TET with LGB.

A Comparison of QM/MM Simulations with and without the Drude Oscillator Model Based on Hydration Free Energies of Simple Solutes

Maintaining a proper balance between specific intermolecular interactions and non-specific solvent interactions is of critical importance in molecular simulations, especially when predicting binding affinities or reaction rates in the condensed phase. The most rigorous metric for characterizing solvent affinity are solvation free energies, which correspond to a transfer from the gas phase into solution. Due to the drastic change of the electrostatic environment during this process, it is also a stringent test of polarization response in the model. Here, we employ both the CHARMM fixed charge and polarizable force fields to predict hydration free energies of twelve simple solutes. The resulting classical ensembles are then reweighted to obtain QM/MM hydration free energies using a variety of QM methods, including MP2, Hartree⁻Fock, density functional methods (BLYP, B3LYP, M06-2X) and semi-empirical methods (OM2 and AM1 ). Our simulations test the compatibility of quantum-mechanical methods with molecular-mechanical water models and solute Lennard⁻Jones parameters. In all cases, the resulting QM/MM hydration free energies were inferior to purely classical results, with the QM/MM Drude force field predictions being only marginally better than the QM/MM fixed charge results. In addition, the QM/MM results for different quantum methods are highly divergent, with almost inverted trends for polarizable and fixed charge water models. While this does not necessarily imply deficiencies in the QM models themselves, it underscores the need to develop consistent and balanced QM/MM interactions. Both the QM and the MM component of a QM/MM simulation have to match, in order to avoid artifacts due to biased solute⁻solvent interactions. Finally, we discuss strategies to improve the convergence and efficiency of multi-scale free energy simulations by automatically adapting the molecular-mechanics force field to the target quantum method.

Molecular dynamics based enhanced sampling of collective variables with very large time steps

Enhanced sampling techniques that target a set of collective variables and that use molecular dynamics as the driving engine have seen widespread application in the computational molecular sciences as a means to explore the free-energy landscapes of complex systems. The use of molecular dynamics as the fundamental driver of the sampling requires the introduction of a time step whose magnitude is limited by the fastest motions in a system. While standard multiple time-stepping methods allow larger time steps to be employed for the slower and computationally more expensive forces, the maximum achievable increase in time step is limited by resonance phenomena, which inextricably couple fast and slow motions. Recently, we introduced deterministic and stochastic resonance-free multiple time step algorithms for molecular dynamics that solve this resonance problem and allow ten- to twenty-fold gains in the large time step compared to standard multiple time step algorithms [P. Minary et al., Phys. Rev. Lett. 93, 150201 (2004); B. Leimkuhler et al., Mol. Phys. 111, 3579-3594 (2013)]. These methods are based on the imposition of isokinetic constraints that couple the physical system to Nosé-Hoover chains or Nosé-Hoover Langevin schemes. In this paper, we show how to adapt these methods for collective variable-based enhanced sampling techniques, specifically adiabatic free-energy dynamics/temperature-accelerated molecular dynamics, unified free-energy dynamics, and by extension, metadynamics, thus allowing simulations employing these methods to employ similarly very large time steps. The combination of resonance-free multiple time step integrators with free-energy-based enhanced sampling significantly improves the efficiency of conformational exploration.

Mapping the Drude polarizable force field onto a multipole and induced dipole model

The induced dipole and the classical Drude oscillator represent two major approaches for the explicit inclusion of electronic polarizability into force field-based molecular modeling and simulations. In this work, we explore the equivalency of these two models by comparing condensed phase properties computed using the Drude force field and a multipole and induced dipole (MPID) model. Presented is an approach to map the electrostatic model optimized in the context of the Drude force field onto the MPID model. Condensed phase simulations on water and 15 small model compounds show that without any reparametrization, the MPID model yields properties similar to the Drude force field with both models yielding satisfactory reproduction of a range of experimental values and quantum mechanical data. Our results illustrate that the Drude oscillator model and the point induced dipole model are different representations of essentially the same physical model. However, results indicate the presence of small differences between the use of atomic multipoles and off-center charge sites. Additionally, results on the use of dispersion particle mesh Ewald further support its utility for treating long-range Lennard Jones dispersion contributions in the context of polarizable force fields. The main motivation in demonstrating the transferability of parameters between the Drude and MPID models is that the more than 15 years of development of the Drude polarizable force field can now be used with MPID formalism without the need for dual-thermostat integrators nor self-consistent iterations. This opens up a wide range of new methodological opportunities for polarizable models.

Assessing many-body contributions to intermolecular interactions of the AMOEBA force field using energy decomposition analysis of electronic structure calculations

In this work, we evaluate the accuracy of the classical AMOEBA model for representing many-body interactions, such as polarization, charge transfer, and Pauli repulsion and dispersion, through comparison against an energy decomposition method based on absolutely localized molecular orbitals (ALMO-EDA) for the water trimer and a variety of ion-water systems. When the 2- and 3-body contributions according to the many-body expansion are analyzed for the ion-water trimer systems examined here, the 3-body contributions to Pauli repulsion and dispersion are found to be negligible under ALMO-EDA, thereby supporting the validity of the pairwise-additive approximation in AMOEBA's 14-7 van der Waals term. However AMOEBA shows imperfect cancellation of errors for the missing effects of charge transfer and incorrectness in the distance dependence for polarization when compared with the corresponding ALMO-EDA terms. We trace the larger 2-body followed by 3-body polarization errors to the Thole damping scheme used in AMOEBA, and although the width parameter in Thole damping can be changed to improve agreement with the ALMO-EDA polarization for points about equilibrium, the correct profile of polarization as a function of intermolecular distance cannot be reproduced. The results suggest that there is a need for re-examining the damping and polarization model used in the AMOEBA force field and provide further insights into the formulations of polarizable force fields in general.

Mesh-free hierarchical clustering methods for fast evaluation of electrostatic interactions of point multipoles

Electrostatic interactions involving point multipoles are being increasingly implemented to achieve higher accuracy in molecular simulations. A major drawback of multipolar electrostatics is the increased computational cost. Here we develop and compare two Cartesian tree algorithms which employ Taylor approximations and hierarchical clustering to speed up the evaluation of point multipole interactions. We present results from applying the algorithms to compute the free space Coulomb potential and forces of different sets of interacting point multipoles with different densities. The methods achieve high accuracy and speedup of more than an order of magnitude over direct sum calculations and scale well in parallel.

The strengths and limitations of effective centroid force models explored by studying isotopic effects in liquid water

The development of effective centroid potentials (ECPs) is explored with both the constrained-centroid and quasi-adiabatic force matching using liquid water as a test system. A trajectory integrated with the ECP is free of statistical noises that would be introduced when the centroid potential is approximated on the fly with a finite number of beads. With the reduced cost of ECP, challenging experimental properties can be studied in the spirit of centroid molecular dynamics. The experimental number density of H2O is 0.38% higher than that of D2O. With the ECP, the H2O number density is predicted to be 0.42% higher, when the dispersion term is not refit. After correction of finite size effects, the diffusion constant of H2O is found to be 21% higher than that of D2O, which is in good agreement with the 29.9% higher diffusivity for H2O observed experimentally. Although the ECP is also able to capture the redshifts of both the OH and OD stretching modes in liquid water, there are a number of properties that a classical simulation with the ECP will not be able to recover. For example, the heat capacities of H2O and D2O are predicted to be almost identical and higher than the experimental values. Such a failure is simply a result of not properly treating quantized vibrational energy levels when the trajectory is propagated with classical mechanics. Several limitations of the ECP based approach without bead population reconstruction are discussed.