Seamless integration of GEM, a density based-force field, for QM/MM simulations via LICHEM, Psi4, and Tinker-HP.

Hybrid quantum mechanics/molecular mechanics (QM/MM) simulations have become an essential tool in computational chemistry, particularly for analyzing complex biological and condensed phase systems. Building on this foundation, our work presents a novel implementation of the Gaussian Electrostatic Model (GEM), a polarizable density-based force field, within the QM/MM framework. This advancement provides seamless integration, enabling efficient and optimized QM/GEM calculations in a single step using the LICHEM Code. We have successfully applied our implementation to water dimers and hexamers, demonstrating the ability to handle water systems with varying numbers of water molecules. Moreover, we have extended the application to describe the double proton transfer of the aspartic acid dimer in a box of water, which highlights the method's proficiency in investigating heterogeneous systems. Our implementation offers the flexibility to perform on-the-fly density fitting or to utilize pre-fitted coefficients to estimate exchange and Coulomb contributions. This flexibility enhances efficiency and accuracy in modeling molecular interactions, especially in systems where polarization effects are significant.

Impact of an Ionic Liquid Solution on Horseradish Peroxidase Activity.

Horseradish peroxidase (HRP) is an enzyme that oxidizes pollutants from wastewater. A previous report indicated that peroxidases can have an enhancement in initial enzymatic activity in an aqueous solution of 0.26 M 1-ethyl-3-methylimidazolium ethyl sulfate ([EMIm][EtSO4]) at neutral pH. However, the atomistic details remain elusive. In the enzymatic landscape of HRP, compound II (Cpd II) plays a key role and involves a histidine (H42) residue. Cpd II exists as oxoferryl (2a) or hydroxoferryl (2b(FeIV)) forms, where 2a is the predominantly observed form in experimental studies. Intriguingly, the ferric 2b(FeIII) form seen in synthetic complexes has not been observed in HRP. Here, we have investigated the structure and dynamics of HRP in pure water and aqueous [EMIm][EtSO4] (0.26 M), as well as the reaction mechanism of 2a to 2b conversion using polarizable molecular dynamics (MD) simulations and quantum mechanics/molecular mechanics (QM/MM) calculations. When HRP is solvated in aq [EMIm][EtSO4], the catalytic water displaces, and H42 directly orients over the ferryl moiety, allowing a direct proton transfer (PT) with a significant energy barrier reduction. Conversely, in neat water, the reaction of 2a to 2b follows the previously reported mechanism. We further investigated the deprotonated form of H42. Analysis of the electric fields at the active site indicates that the aq [EMIm][EtSO4] medium facilitates the reaction by providing a more favorable environment compared with the system solvated in neat water. Overall, the atomic level supports the previous experimental observations and underscores the importance of favorable electric fields in the active site to promote catalysis.

Relative cooperativity in neutral and charged molecular clusters using QM/MM calculations.

QM/MM methods have been used to study electronic structure properties and chemical reactivity in complex molecular systems where direct electronic structure calculations are not feasible. In our previous work, we showed that non-polarizable force fields, by design, describe intermolecular interactions through pairwise interactions, overlooking many-body interactions involving three or more particles. In contrast, polarizable force fields account partially for many-body effects through polarization, but still handle van der Waals and permanent electrostatic interactions pairwise. We showed that despite those limitations, polarizable and non-polarizable force fields can reproduce relative cooperativity achieved using density functional theory due to error compensation mechanisms. In this contribution, we assess the performance of QM/MM methods in reproducing these phenomena. Our study highlights the significance of the QM region size and force field choice in QM/MM calculations, emphasizing the importance of parameter validation to obtain accurate interaction energy predictions.

Cooperativity and Frustration Effects (or Lack Thereof) in Polarizable and Non-polarizable Force Fields.

Understanding cooperativity and frustration is crucial for studying biological processes such as molecular recognition and protein aggregation. Force fields have been extensively utilized to explore cooperativity in the formation of protein secondary structures and self-assembled systems. Multiple studies have demonstrated that polarizable force fields provide more accurate descriptions of this phenomenon compared to fixed-charge pairwise nonpolarizable force fields, thanks to the incorporation of polarization effects. In this study, we assess the performance of the AMOEBA polarizable force field and the AMBER and OPLS nonpolarizable pairwise force fields in capturing positive and negative cooperativity recently explored in neutral and charged molecular clusters using density functional theory. Our findings show that polarizable and nonpolarizable force fields qualitatively reproduce the relative cooperativity observed in electron structure calculations. However, AMBER and OPLS fail to describe absolute cooperativity. In contrast, AMOEBA accounts for the absolute cooperativity by considering interactions beyond pairwise interactions. According to the energy decomposition analysis, it is observed that the electrostatic interactions calculated with the AMBER and OPLS force fields seem to play an important and counterintuitive role in reproducing the adiabatic interaction energies calculated with density functional theory. However, it is important to note that these force fields, due to their nature, do not explicitly incorporate many-body effects, which limits their ability to accurately describe cooperativity. On the other hand, frustration in polarizable and nonpolarizable force fields is caused by changes in bond stretching and angle bending terms of the building blocks when they are forming a complex.

Polarizable MD and QM/MM investigation of acrylamide-based leads to target the main protease of SARS-CoV-2.

The main protease (Mpro) of SARS-CoV-2 is an essential enzyme for the replication of the virus causing the COVID-19 pandemic. Because there is no known homologue in humans, it has been proposed as a primary target for antiviral drug development. Here, we explore the potential of five acrylamide-based molecules as possible covalent inhibitors, leading to target MPro by docking, followed by polarizable molecular dynamics (MD) and quantum mechanics/molecular mechanics (QM/MM) calculations. All calculations involving a classical potential were calculated with the AMOEBABIO18 polarizable force field, while electronic structure calculations were performed within the framework of density functional theory. Selected docking poses for each of the five compounds were used for MD simulations, which suggest only one of the tested leads remains bound in a catalytically active orientation. The QM/MM results for the covalent attachment of the promising lead to the catalytic serine suggest that this process is thermodynamically feasible but kinetically unlikely. Overall, our results are consistent with the low labeling percentages determined experimentally and may be useful for further development of acrylamide-based leads.

Electronic Structure Effects Related to the Origin of the Remarkable Near-Infrared Absorption of Blastochloris viridis' Light Harvesting 1-Reaction Center Complex.

Various photosynthetic organisms have evolved to absorb light in different regions of the visible light spectrum, thus adapting to the various lighting conditions available on Earth. While most of these autotrophic organisms absorb wavelengths around the 700-800 nm region, some are capable of red-shifted absorptions above this range, but none as remarkably as Blastochloris viridis whose main absorption is observed at 1015 nm, approximately 220 nm (0.34 eV) lower in energy than their main constituent pigments, BChl-b, whose main absorption is observed at 795 nm. The structure of its light harvesting 1-reaction center was recently elucidated by cryo-EM; however, the electronic structure details behind this red-shifted absorption remain unattended. We used hybrid quantum mechanics/molecular mechanics (QM/MM) calculations to optimize one of the active centers and performed classical molecular dynamics (MD) simulations to sample conformations beyond the optimized structure. We did excited state calculations with the time-dependent density functional theory method at the CAM-B3LYP/cc-pVDZ level of theory. We reproduced the near IR absorption by sequentially modifying the number of components involved in our systems using representative structures from the calculated MD ensemble. Natural transition orbital analysis reveals the participation of the BChl-b fragments to the main transition in the native structure and the structures obtained from the QM/MM and MD simulations. H-bonding pigment-protein interactions play a role on the conformation stabilization and orientation; however, the bacteriochlorin ring conformations and the exciton delocalization are the most relevant factors to explain the red-shifting phenomenon.

Structural and electronic analysis of the octarepeat region of prion protein with four Cu2+ by polarizable MD and QM/MM simulations.

Prions have been linked to neurodegenerative diseases that affect various species of mammals including humans. The prion protein, located mainly in neurons, is believed to play the role of metal ion transporter. High levels of copper ions have been related to structural changes. A 32-residue region of the N-terminal domain, known as octarepeat, can bind up to four copper ions. Different coordination modes have been observed and are strongly dependent on Cu2+ concentration. Many theoretical studies carried out so far have focused on studying the coordination modes of a single copper ion. In this work we investigate the octarepeat region coordinated with four copper ions. Molecular dynamics (MD) and hybrid quantum mechanics/molecular mechanics (QM/MM) simulations using the polarizable AMOEBA force field have been carried out. The polarizable MD simulations starting from a fully extended conformation indicate that the tetra-Cu2+/octarepeat complex forms a globular structure. The globular form is stabilized by interactions between Cu2+ and tryptophan residues resulting in some coordination sites observed to be in close proximity, in agreement with experimental results. Subsequent QM/MM simulations on several snapshots suggests the system is in a high-spin quintet state, with all Cu2+ bearing one single electron, and all unpaired electrons are ferromagnetically coupled. NMR simulations on selected structures provides insights on the chemical shifts of the first shell ligands around the metals with respect to inter-metal distances.

Development and application of quantum mechanics/molecular mechanics methods with advanced polarizable potentials.

Quantum mechanics/molecular mechanics (QM/MM) simulations are a popular approach to study various features of large systems. A common application of QM/MM calculations is in the investigation of reaction mechanisms in condensed-phase and biological systems. The combination of QM and MM methods to represent a system gives rise to several challenges that need to be addressed. The increase in computational speed has allowed the expanded use of more complicated and accurate methods for both QM and MM simulations. Here, we review some approaches that address several common challenges encountered in QM/MM simulations with advanced polarizable potentials, from methods to account for boundary across covalent bonds and long-range effects, to polarization and advanced embedding potentials.

Origin of cooperativity in hydrogen bonding.

The origin of non-additivity in hydrogen bonds (H-bonds), usually termed as H-bond cooperativity, is investigated in H-bonded linear chains. It is shown that H-bond cooperativity originates solely from classical electrostatics. The latter is corroborated by comparing the H-bond cooperativity in infinitely-long H-bonded hydrogen cyanide, 4-pyridone and formamide chains, assessed using density functional theory (DFT), against the strengthening of the dipole-dipole interaction upon the formation of an infinite chain of effective point-dipoles. It is found that the magnitude of these effective point-dipoles is a consequence of mutual polarization and additional effects beyond a polarizable point-dipole model. Nevertheless, the effective point-dipoles are fully determined once a single H-bond is formed, indicating that quantum effects involved in H-bonding are circumscribed to nearest-neighbor interactions only; i.e. in a linear chain of H-bonds, quantum effects do not contribute to the H-bond non-additivity. This finding is verified by estimating cooperativity along the dissociation path of H-bonds in the infinite chains, using two empirical parameters that account for polarizability, together with DFT association energies and molecular dipoles of solely monomers and dimers.

On cooperative effects and aggregation of GNNQQNY and NNQQNY peptides.

Some health disturbances like neurodegenerative diseases are associated to the presence of amyloids. GNNQQNY and NNQQNY peptides are considered as prototypical examples for studying the formation of amyloids. These exhibit quite different aggregation behaviors despite they solely differ in size by one residue. To get insight into the reasons for such difference, we have examined association energies of aggregates (parallel ß-sheets, fibril-spines, and crystal structures) from GNNQQNY and NNQQY using density functional theory. As we found that GNNQQNY tends to form a zwitterion in the crystal structure, we have investigated the energetics of parallel ß-sheets and fibril-spines in the canonical and zwitterionic states. We found that the formation of GNNQQNY aggregates is energetically more favored than the formation of the NNQQNY ones. We show that the latter is connected to the network of hydrogen bonds formed by each aggregate. Moreover, we found that the formation of some NNQQNY aggregates is anticooperative, whereas cooperative with GNNQQNY. These results have interesting implications for deciphering the factors determining peptide aggregation propensities.