Enhanced semiempirical QM methods for biomolecular interactions

dxtb-An efficient and fully differentiable framework for extended tight-binding

Automatic differentiation (AD) emerged as an integral part of machine learning, accelerating model development by enabling gradient-based optimization without explicit analytical derivatives. Recently, the benefits of AD and computing arbitrary-order derivatives with respect to any variable were also recognized in the field of quantum chemistry. In this work, we present dxtb-an open-source, fully differentiable framework for semiempirical extended tight-binding (xTB) methods. Developed entirely in Python and leveraging PyTorch for array operations, dxtb facilitates extensibility and rapid prototyping while maintaining computational efficiency. Through comprehensive code vectorization and optimization, we essentially reach the speed of compiled xTB programs for high-throughput calculations of small molecules. The excellent performance also scales to large systems, and batch operability yields additional benefits for execution on parallel hardware. In particular, energy evaluations are on par with existing programs, whereas the speed of automatically differentiated nuclear derivatives is only 2 to 5 times slower compared to their analytical counterparts. We showcase the utility of AD in dxtb by calculating various molecular and spectroscopic properties, highlighting its capacity to enhance and simplify such evaluations. Furthermore, the framework streamlines optimization tasks and offers seamless integration of semiempirical quantum chemistry in machine learning, paving the way for physics-inspired end-to-end differentiable models. Ultimately, dxtb aims to further advance the capabilities of semiempirical methods, providing an extensible foundation for future developments and hybrid machine learning applications. The framework is accessible at https://github.com/grimme-lab/dxtb.

Development of NOTCH, an all-electron, beyond-NDDO semiempirical method: Application to diatomic molecules

In this work, we develop a new semiempirical method, dubbed NOTCH (Natural Orbital Tied Constructed Hamiltonian). Compared to existing semiempirical methods, NOTCH is less empirical in its functional form as well as parameterization. Specifically, in NOTCH, (1) the core electrons are treated explicitly; (2) the nuclear-nuclear repulsion term is calculated analytically, without any empirical parameterization; (3) the contraction coefficients of the atomic orbital (AO) basis depend on the coordinates of the neighboring atoms, which allows the size of AOs to depend on the molecular environment, despite the fact that a minimal basis set is used; (4) the one-center integrals of free atoms are derived from scalar relativistic multireference equation-of-motion coupled cluster calculations instead of empirical fitting, drastically reducing the number of necessary empirical parameters; (5) the (AA|AB) and (AB|AB)-type two-center integrals are explicitly included, going beyond the neglect of differential diatomic overlap approximation; and (6) the integrals depend on the atomic charges, effectively mimicking the "breathing" of AOs when the atomic charge varies. For this preliminary report, the model has been parameterized for the elements H-Ne, giving only 8 empirical global parameters. Preliminary results on the ionization potentials, electron affinities, and excitation energies of atoms and diatomic molecules, as well as the equilibrium geometries, vibrational frequencies dipole moments, and bond dissociation energies of diatomic molecules, show that the accuracy of NOTCH rivals or exceeds those of popular semiempirical methods (including PM3, PM7, OM2, OM3, GFN-xTB, and GFN2-xTB) as well as the cost-effective ab initio method Hartree-Fock-3c.

ReaxFF-based nonadiabatic dynamics method for azobenzene derivatives

ReaxFF reactive force fields have been parameterized for the ground and first excited states of azobenzene and its derivatives. In addition, an extended set of ab initio reference data ensures wide applicability, including to azosystems in complex environments. Based on the optimized force fields, nonadiabatic surface hopping simulations produce photoisomerization quantum yields and decay times of azobenzene, both in the gas phase and in n-hexane solution, in reasonable agreement with higher level theory and experiment. The transferability to other azo-compounds is illustrated for different arylazopyrazoles as well as ethylene-bridged azobenzene. Moreover, it has been shown that the model can be easily extended to adsorbates on metal surfaces. The simulation of the ring-opening of cyclobutene triggered by the photoisomerization of azobenzene in a macrocycle highlights the advantages of a reactive force field model.

Effect of the Solvate Environment of Lithium Cations on the Resistance of the Polymer Electrolyte/Electrode Interface in a Solid-State Lithium Battery

The effect of the composition of liquid electrolytes in the bulk and at the interface with the LiFePO₄ cathode on the operation of a solid-state lithium battery with a nanocomposite polymer gel electrolyte based on polyethylene glycol diacrylate and SiO₂ was studied. The self-diffusion coefficients on the 7Li, 1H, and 19F nuclei in electrolytes based on LiBF₄ and LiTFSI salts in solvents (gamma-butyrolactone, dioxolane, dimethoxyethane) were measured by nuclear magnetic resonance (NMR) with a magnetic field gradient. Four compositions of the complex electrolyte system were studied by high-resolution NMR. The experimentally obtained ¹H chemical shifts are compared with those theoretically calculated by quantum chemical modeling. This made it possible to suggest the solvate shell compositions that facilitate the rapid transfer of the Li⁺ cation at the nanocomposite electrolyte/LiFePO₄ interface and ensure the stable operation of a solid-state lithium battery.

ULYSSES: An Efficient and Easy to Use Semiempirical Library for C+

We introduce ULYSSES, a user-friendly and robust C++ library for semiempirical quantum chemical calculations. In its current version, ULYSSES is equipped with a large set of different semiempirical models, most of which are based on the Neglect of Diatomic Differential Overlap (NDDO) approximation. Empirical corrections for dispersion and hydrogen bonding are available for most methods, so that higher quality is achieved in the calculation of energies of nonbonded complexes. The library is furthermore equipped with geometry optimization, as well as modules for calculating molecular properties of general interest. Ideal gas thermodynamics is available and allows single structure as well as conformer (multistructure) averaged properties to be calculated. We offer the possibility to use several vibrational partition functions according to the nature of interactions being studied: for covalent systems, the traditional harmonic oscillator approximation is available; for nonbonded complexes, we systematically extended the partition function proposed by Grimme for all thermodynamic functions. The library is also capable of running Born-Oppenheimer molecular dynamics.

Fast and Accurate Quantum Mechanical Modeling of Large Molecular Systems Using Small Basis Set Hartree-Fock Methods Corrected with Atom-Centered Potentials

There has been significant interest in developing fast and accurate quantum mechanical methods for modeling large molecular systems. In this work, by utilizing a machine learning regression technique, we have developed new low-cost quantum mechanical approaches to model large molecular systems. The developed approaches rely on using one-electron Gaussian-type functions called atom-centered potentials (ACPs) to correct for the basis set incompleteness and the lack of correlation effects in the underlying minimal or small basis set Hartree-Fock (HF) methods. In particular, ACPs are proposed for ten elements common in organic and bioorganic chemistry (H, B, C, N, O, F, Si, P, S, and Cl) and four different base methods: two minimal basis sets (MINIs and MINIX) plus a double-ζ basis set (6-31G*) in combination with dispersion-corrected HF (HF-D3/MINIs, HF-D3/MINIX, HF-D3/6-31G*) and the HF-3c method. The new ACPs are trained on a very large set (73 832 data points) of noncovalent properties (interaction and conformational energies) and validated additionally on a set of 32 048 data points. All reference data are of complete basis set coupled-cluster quality, mostly CCSD(T)/CBS. The proposed ACP-corrected methods are shown to give errors in the tenths of a kcal/mol range for noncovalent interaction energies and up to 2 kcal/mol for molecular conformational energies. More importantly, the average errors are similar in the training and validation sets, confirming the robustness and applicability of these methods outside the boundaries of the training set. In addition, the performance of the new ACP-corrected methods is similar to complete basis set density functional theory (DFT) but at a cost that is orders of magnitude lower, and the proposed ACPs can be used in any computational chemistry program that supports effective-core potentials without modification. It is also shown that ACPs improve the description of covalent and noncovalent bond geometries of the underlying methods and that the improvement brought about by the application of the ACPs is directly related to the number of atoms to which they are applied, allowing the treatment of systems containing some atoms for which ACPs are not available. Overall, the ACP-corrected methods proposed in this work constitute an alternative accurate, economical, and reliable quantum mechanical approach to describe the geometries, interaction energies, and conformational energies of systems with hundreds to thousands of atoms.

Thermochemical and Quantum Descriptor Calculations for Gaining Insight into Ricin Toxin A (RTA) Inhibitors

In this work, we performed a study to assess the interactions between the ricin toxin A (RTA) subunit of ricin and some of its inhibitors using modern semiempirical quantum chemistry and ONIOM quantum mechanics/molecular mechanics (QM/MM) methods. Two approaches were followed (calculation of binding enthalpies, ΔH _bind, and reactivity quantum chemical descriptors) and compared with the respective half-maximal inhibitory concentration (IC₅₀) experimental data, to gain insight into RTA inhibitors and verify which quantum chemical method would better describe RTA-ligand interactions. The geometries for all RTA-ligand complexes were obtained after running classical molecular dynamics simulations in aqueous media. We found that single-point energy calculations of ΔH _bind with the PM6-DH+, PM6-D3H4, and PM7 semiempirical methods and ONIOM QM/MM presented a good correlation with the IC₅₀ data. We also observed, however, that the correlation decreased significantly when we calculated ΔH _bind after full-atom geometry optimization with all semiempirical methods. Based on the results from reactivity descriptors calculations for the cases studied, we noted that both types of interactions, molecular overlap and electrostatic interactions, play significant roles in the overall affinity of these ligands for the RTA binding pocket.

Antioxidant and antimicrobial material by grafting of L-arginine onto enzymatic poly(gallic acid)

Microwave-mediated grafting of L-Arg onto naturally derived and stable multiradical poly(gallic acid) (PGAL) in aqueous media has been successfully achieved. This polymeric material has no adverse effect in human cells as there is no hemolytic activity upon MTT and Neutral Red assays. The analytical and computational characterization studies carried out in this study describe a helical molecular structure with random incorporation of L-Arginine pendant groups from PGAL's backbone. The antioxidant properties of the precursor polymer are preserved as proved by the elimination of stable DPPH and hydroxyl radical scavenging, as well as the FRAP and ORAC assays. Regarding the latter, the oxygen radical inhibition is enhanced compared to PGAL, which is attributed to the guanidyl moieties. PGAL-g-L-Arg displays antimicrobial activity against Gram (+) Listeria monocytogenes and Staphylococcus aureus strains with a MIC of 0.8 g/L and a bacteriostatic effect against Gram (-) Escherichia coli. Additionally, scanning electron and confocal fluorescence microscopies as well as crystal violet colorimetric assay demonstrate that the mechanism involved in the bacterial inhibition is related to the formation of porous channels on the membrane, which is discussed according to the helical secondary structure of the polymer and the amino acid guanidyl moieties interacting to bacterial membranes.

Robust Atomistic Modeling of Materials, Organometallic, and Biochemical Systems

Modern chemistry seems to be unlimited in molecular size and elemental composition. Metal‐organic frameworks or biological macromolecules involve complex architectures and a large variety of elements. Yet, a general and broadly applicable theoretical method to describe the structures and interactions of molecules beyond the 1000‐atom size regime semi‐quantitatively is not self‐evident. For this purpose, a generic force field named GFN‐FF is presented, which is completely newly developed to enable fast structure optimizations and molecular‐dynamics simulations for basically any chemical structure consisting of elements up to radon. The freely available computer program requires only starting coordinates and elemental composition as input from which, fully automatically, all potential‐energy terms are constructed. GFN‐FF outperforms other force fields in terms of generality and accuracy, approaching the performance of much more elaborate quantum‐mechanical methods in many cases.

Revisiting the Potential Energy Surface of the Stacked Cytosine Dimer: FNO-CCSD(T) Interaction Energies, SAPT Decompositions, and Benchmarking

Nucleobase stacking interactions are crucial for the stability of nucleic acids. This study investigates base stacking energies of the cytosine homodimer in different configurations, including intermolecular separation plots, detailed twist dependence, and displaced structures. Highly accurate ab initio quantum chemical single point energies using an energy function based on MP2 complete basis set extrapolation ([6 → 7]ZaPa-NR) and a CCSD(T)/cc-pVTZ-F12 high-level correction are presented as new reference data, providing the most accurate stacking energies of nucleobase dimers currently available. Accurate SAPT2+(3)δMP2 energy decomposition is used to obtain detailed insights into the nature of base stacking interactions at varying vertical distances and twist values. The ab initio symmetry adapted perturbation theory (SAPT) energy decomposition suggests that the base stacking originates from an intricate interplay between dispersion attraction, short-range exchange-repulsion, and Coulomb interaction. The interpretation of the SAPT data is a complex issue as key energy terms vary substantially in the region of optimal (low energy) base stacking geometries. Thus, attempts to highlight one leading stabilizing SAPT base stacking term may be misleading and the outcome strongly depends on the used geometries within the range of geometries sampled in nucleic acids upon thermal fluctuations. Modern dispersion-corrected density functional theory (among them DSD-BLYP-D3, ωB97M-V, and ωB97M-D3BJ) is benchmarked and often reaches up to spectroscopic accuracy (below 1 kJ/mol). The classical AMBER force field is benchmarked with multiple different sets of point-charges (e.g. HF, DFT, and MP2-based) and is found to produce reasonable agreement with the benchmark data.

In silico identification of genetic mutations conferring resistance to acetohydroxyacid synthase inhibitors: A case study of Kochia scoparia

Mutations that confer herbicide resistance are a primary concern for herbicide-based chemical control of invasive plants and are often under-characterized structurally and functionally. As the outcome of selection pressure, resistance mutations usually result from repeated long-term applications of herbicides with the same mode of action and are discovered through extensive field trials. Here we used acetohydroxyacid synthase (AHAS) of Kochia scoparia (KsAHAS) as an example to demonstrate that, given the sequence of a target protein, the impact of genetic mutations on ligand binding could be evaluated and resistance mutations could be identified using a biophysics-based computational approach. Briefly, the 3D structures of wild-type (WT) and mutated KsAHAS-herbicide complexes were constructed by homology modeling, docking and molecular dynamics simulation. The resistance profile of two AHAS-inhibiting herbicides, tribenuron methyl and thifensulfuron methyl, was obtained by estimating their binding affinity with 29 KsAHAS (1 WT and 28 mutated) using 6 molecular mechanical (MM) and 18 hybrid quantum mechanical/molecular mechanical (QM/MM) methods in combination with three structure sampling strategies. By comparing predicted resistance with experimentally determined resistance in the 29 biotypes of K. scoparia field populations, we identified the best method (i.e., MM-PBSA with single structure) out of all tested methods for the herbicide-KsAHAS system, which exhibited the highest accuracy (up to 100%) in discerning mutations conferring resistance or susceptibility to the two AHAS inhibitors. Our results suggest that the in silico approach has the potential to be widely adopted for assessing mutation-endowed herbicide resistance on a case-by-case basis.

Semiempirical Quantum-Chemical Methods with Orthogonalization and Dispersion Corrections

We present two new semiempirical quantum-chemical methods with orthogonalization and dispersion corrections: ODM2 and ODM3 (ODM x). They employ the same electronic structure model as the OM2 and OM3 (OM x) methods, respectively. In addition, they include Grimme's dispersion correction D3 with Becke-Johnson damping and three-body corrections E ABC for Axilrod-Teller-Muto dispersion interactions as integral parts. Heats of formation are determined by adding explicitly computed zero-point vibrational energy and thermal corrections, in contrast to standard MNDO-type and OM x methods. We report ODM x parameters for hydrogen, carbon, nitrogen, oxygen, and fluorine that are optimized with regard to a wide range of carefully chosen state-of-the-art reference data. Extensive benchmarks show that the ODM x methods generally perform better than the available MNDO-type and OM x methods for ground-state and excited-state properties, while they describe noncovalent interactions with similar accuracy as OM x methods with a posteriori dispersion corrections.

Improving the accuracy of predicting protein-ligand binding-free energy with semiempirical quantum chemistry charge

Aim: It is a challenge to predict binding-free energy (ΔG) accurately. Methodology/results: For accurate ΔG prediction, a new strategy combining solvated interaction energy (SIE) or molecular mechanics/generalized Born surface area (MM/GBSA) approach with the Coulson charge of both protein and ligand calculated by semiempirical quantum mechanics (SQM), named SIE-SQMPC or MM/GBSA-SQMPC approach, was developed and tested on 50 protein-ligand complexes. Both approaches achieved higher correlation (R ²) between experimental and predicted ΔG than that with Amber-ff03 charge, even for ligands with highly different scaffolds. But, SIE-SQMPC is computationally much faster than MM/GBSA-SQMPC. Conclusion: SIE-SQMPC provided an effective alternative to predict ΔG of protein-ligand binding (R ² = 0.66-0.94 for SIE-AM1; R ² = 0.59-0.98 for SIE-PM7), which has the potential of high-throughput processing for molecular docking and drug design.

GFN2-xTB-An Accurate and Broadly Parametrized Self-Consistent Tight-Binding Quantum Chemical Method with Multipole Electrostatics and Density-Dependent Dispersion Contributions

An extended semiempirical tight-binding model is presented, which is primarily designed for the fast calculation of structures and noncovalent interaction energies for molecular systems with roughly 1000 atoms. The essential novelty in this so-called GFN2-xTB method is the inclusion of anisotropic second order density fluctuation effects via short-range damped interactions of cumulative atomic multipole moments. Without noticeable increase in the computational demands, this results in a less empirical and overall more physically sound method, which does not require any classical halogen or hydrogen bonding corrections and which relies solely on global and element-specific parameters (available up to radon, Z = 86). Moreover, the atomic partial charge dependent D4 London dispersion model is incorporated self-consistently, which can be naturally obtained in a tight-binding picture from second order density fluctuations. Fully analytical and numerically precise gradients (nuclear forces) are implemented. The accuracy of the method is benchmarked for a wide variety of systems and compared with other semiempirical methods. Along with excellent performance for the "target" properties, we also find lower errors for "off-target" properties such as barrier heights and molecular dipole moments. High computational efficiency along with the improved physics compared to its precursor GFN-xTB makes this method well-suited to explore the conformational space of molecular systems. Significant improvements are furthermore observed for various benchmark sets, which are prototypical for biomolecular systems in aqueous solution.

Affinity Calculations of Cyclodextrin Host-Guest Complexes: Assessment of Strengths and Weaknesses of End-Point Free Energy Methods

The end-point methods like MM/PBSA or MM/GBSA estimate the free energy of a biomolecule by combining its molecular mechanics energy with solvation free energy and entropy terms. On the one hand, their performance largely depends on the particular system of interest, and despite numerous attempts to improve their reliability that have resulted in many variants, there is still no clear alternative to improve their accuracy. On the other hand, the relatively small cyclodextrin host-guest complexes, for which high-quality binding calorimetric data are usually available, are becoming reference models for testing the accuracy of free energy methods. In this work, we further assess the performance of various MM/PBSA-like approaches as applied to cyclodextrin complexes. To this end, we select a set of complexes between β-cyclodextrin and 57 small organic molecules that has been previously studied with the binding energy distribution analysis method in combination with an implicit solvent model ( Wickstrom, L.; He, P.; Gallicchio, E.; Levy, R. M. J. Chem. Theory Comput. 2013 , 9 , 3136 - 3150 ). For each complex, a conventional 1.0 μs molecular dynamics simulation in explicit solvent is performed. Then we employ semiempirical quantum chemical calculations, several variants of the MM-PB(GB)SA methods, entropy estimations, etc., to assess the reliability of the end-point affinity calculations. The best end-point protocol in this study, which combines DFTB3 energies with entropy corrections, yields estimations of the binding free energies that still have substantial errors (RMSE = 2.2 kcal/mol), but it exhibits a good prediction capacity in terms of ligand ranking ( R² = 0.66) that is close to or even better than that of rigorous free energy methodologies. Our results can be helpful to discriminate between the intrinsic limitations of the end-point methods and other sources of error, such as the underlying energy and continuum solvation methods.

QM/MM Calculations on Protein-RNA Complexes: Understanding Limitations of Classical MD Simulations and Search for Reliable Cost-Effective QM Methods

Although atomistic explicit-solvent Molecular Dynamics (MD) is a popular tool to study protein-RNA recognition, satisfactory MD description of protein-RNA complexes is not always achieved. Unfortunately, it is often difficult to separate MD simulation instabilities primarily caused by the simple point-charge molecular mechanics (MM) force fields from problems related to the notorious uncertainties in the starting structures. Herein, we report a series of large-scale QM/MM calculations on the U1A protein-RNA complex. This experimentally well-characterized system has an intricate protein-RNA interface, which is very unstable in MD simulations. The QM/MM calculations identify several H-bonds poorly described by the MM method and thus indicate the sources of instabilities of the U1A interface in MD simulations. The results suggest that advanced QM/MM computations could be used to indirectly rationalize problems seen in MM-based MD simulations of protein-RNA complexes. As the most accurate QM method, we employ the computationally demanding meta-GGA density functional TPSS-D3(BJ)/def2-TZVP level of theory. Because considerably faster methods would be needed to extend sampling and to study even larger protein-RNA interfaces, a set of low-cost QM/MM methods is compared to the TPSS-D3(BJ)/def2-TZVP data. The PBEh-3c and B97-3c density functional composite methods appear to be suitable for protein-RNA interfaces. In contrast, HF-3c and the tight-binding Hamiltonians DFTB3-D3 and GFN-xTB perform unsatisfactorily and do not provide any advantage over the MM description. These conclusions are supported also by similar analysis of a simple HutP protein-RNA interface, which is well-described by MD with the exception of just one H-bond. Some other methodological aspects of QM/MM calculations on protein-RNA interfaces are discussed.

Computational Biochemistry-Enzyme Mechanisms Explored

Understanding enzyme mechanisms is a major task to achieve in order to comprehend how living cells work. Recent advances in biomolecular research provide huge amount of data on enzyme kinetics and structure. The analysis of diverse experimental results and their combination into an overall picture is, however, often challenging. Microscopic details of the enzymatic processes are often anticipated based on several hints from macroscopic experimental data. Computational biochemistry aims at creation of a computational model of an enzyme in order to explain microscopic details of the catalytic process and reproduce or predict macroscopic experimental findings. Results of such computations are in part complementary to experimental data and provide an explanation of a biochemical process at the microscopic level. In order to evaluate the mechanism of an enzyme, a structural model is constructed which can be analyzed by several theoretical approaches. Several simulation methods can and should be combined to get a reliable picture of the process of interest. Furthermore, abstract models of biological systems can be constructed combining computational and experimental data. In this review, we discuss structural computational models of enzymatic systems. We first discuss various models to simulate enzyme catalysis. Furthermore, we review various approaches how to characterize the enzyme mechanism both qualitatively and quantitatively using different modeling approaches.

Approximate quantum chemical methods for modelling carbohydrate conformation and aromatic interactions: β-cyclodextrin and its adsorption on a single-layer graphene sheet

Adsorption of carbohydrates on graphene has the potential to improve graphene dispersibility in water. Here we assess the ability of DFTB-based and NDDO-based quantum chemical methods to model β-cyclodextrin conformations and interactions with graphene.

Explicit consideration of spatial hydrogen bonding direction for activity coefficient prediction based on implicit solvation calculations

Directional hydrogen bonding is introduced to implicit solvation calculations for improved prediction of solvation properties and phase equilibria of associating fluids.

Recent Progress in Treating Protein-Ligand Interactions with Quantum-Mechanical Methods

We review the first successes and failures of a “new wave” of quantum chemistry-based approaches to the treatment of protein/ligand interactions. These approaches share the use of “enhanced”, dispersion (D), and/or hydrogen-bond (H) corrected density functional theory (DFT) or semi-empirical quantum mechanical (SQM) methods, in combination with ensemble weighting techniques of some form to capture entropic effects. Benchmark and model system calculations in comparison to high-level theoretical as well as experimental references have shown that both DFT-D (dispersion-corrected density functional theory) and SQM-DH (dispersion and hydrogen bond-corrected semi-empirical quantum mechanical) perform much more accurately than older DFT and SQM approaches and also standard docking methods. In addition, DFT-D might soon become and SQM-DH already is fast enough to compute a large number of binding modes of comparably large protein/ligand complexes, thus allowing for a more accurate assessment of entropic effects.

Semiempirical Quantum Mechanical Methods for Noncovalent Interactions for Chemical and Biochemical Applications

Semiempirical (SE) methods can be derived from either Hartree-Fock or density functional theory by applying systematic approximations, leading to efficient computational schemes that are several orders of magnitude faster than ab initio calculations. Such numerical efficiency, in combination with modern computational facilities and linear scaling algorithms, allows application of SE methods to very large molecular systems with extensive conformational sampling. To reliably model the structure, dynamics, and reactivity of biological and other soft matter systems, however, good accuracy for the description of noncovalent interactions is required. In this review, we analyze popular SE approaches in terms of their ability to model noncovalent interactions, especially in the context of describing biomolecules, water solution, and organic materials. We discuss the most significant errors and proposed correction schemes, and we review their performance using standard test sets of molecular systems for quantum chemical methods and several recent applications. The general goal is to highlight both the value and limitations of SE methods and stimulate further developments that allow them to effectively complement ab initio methods in the analysis of complex molecular systems.

Determination of the protonation state for the catalytic dyad in β-secretase when bound to hydroxyethylamine transition state analogue inhibitors: A molecular dynamics simulation study

BACE1 is an aspartyl protease of pharmacological interest for its direct participation in Alzheimer's disease (AD) through β-amyloid peptide production. Two aspartic acid residues are present in the BACE1 catalytic region which can adopt multiple protonation states depending on the chemical nature of its inhibitors, i.e., monoprotonated, diprotonated and di-deprotonated states. In the present study a series of protein-ligand molecular dynamics (MD) simulations was carried out to identify the most feasible protonation state adopted by the catalytic dyad in the presence of hydroxyethylamine transition state analogue inhibitors. The MD trajectories revealed that the di-deprotonated state is most prefered in the presence of hydroxyethilamine (HEA) family inhibitors. This appears as a result after evaluating, for all 9 protonation state configurations during the simulation time, the deviations of a set of distances and dihedral angles measured on the ligand, protein and protein-ligand complex with reference to an X-ray experimental BACE1/HEA crystallographic structure. These results will help to clarify the phenomena related to the HEAs inhibitory pathway, and improve HEAs databases' virtual screening and ligand design processes targeting β-secretase protein.

Semiempirical Quantum-Chemical Orthogonalization-Corrected Methods: Theory, Implementation, and Parameters

Semiempirical orthogonalization-corrected methods (OM1, OM2, and OM3) go beyond the standard MNDO model by explicitly including additional interactions into the Fock matrix in an approximate manner (Pauli repulsion, penetration effects, and core–valence interactions), which yields systematic improvements both for ground-state and excited-state properties. In this Article, we describe the underlying theoretical formalism of the OMx methods and their implementation in full detail, and we report all relevant OMx parameters for hydrogen, carbon, nitrogen, oxygen, and fluorine. For a standard set of mostly organic molecules commonly used in semiempirical method development, the OMx results are found to be superior to those from standard MNDO-type methods. Parametrized Grimme-type dispersion corrections can be added to OM2 and OM3 energies to provide a realistic treatment of noncovalent interaction energies, as demonstrated for the complexes in the S22 and S66×8 test sets.

Predicting the relative binding affinity of mineralocorticoid receptor antagonists by density functional methods

In drug discovery, prediction of binding affinity ahead of synthesis to aid compound prioritization is still hampered by the low throughput of the more accurate methods and the lack of general pertinence of one method that fits all systems. Here we show the applicability of a method based on density functional theory using core fragments and a protein model with only the first shell residues surrounding the core, to predict relative binding affinity of a matched series of mineralocorticoid receptor (MR) antagonists. Antagonists of MR are used for treatment of chronic heart failure and hypertension. Marketed MR antagonists, spironolactone and eplerenone, are also believed to be highly efficacious in treatment of chronic kidney disease in diabetes patients, but is contra-indicated due to the increased risk for hyperkalemia. These findings and a significant unmet medical need among patients with chronic kidney disease continues to stimulate efforts in the discovery of new MR antagonist with maintained efficacy but low or no risk for hyperkalemia. Applied on a matched series of MR antagonists the quantum mechanical based method gave an R(2) = 0.76 for the experimental lipophilic ligand efficiency versus relative predicted binding affinity calculated with the M06-2X functional in gas phase and an R(2) = 0.64 for experimental binding affinity versus relative predicted binding affinity calculated with the M06-2X functional including an implicit solvation model. The quantum mechanical approach using core fragments was compared to free energy perturbation calculations using the full sized compound structures.

A computational study of the self-assembly of the RFFFR peptide

The β-amyloid peptide sequence, LVFFA, inspired the investigation of the fiber formation potential of the RFFFR peptide. The self-assembly was studied in silico by coarse grained-, atomistic molecular dynamics simulations and semi-empirical quantum mechanical calculations. The fiber formation was found to occur according to a three step process starting with the emergence of small aggregates that join together and form fiber segments that eventually form one continuous fiber. From a series of simulations the critical fiber concentration was determined to be in the interval between 70 mM and 100 mM. To obtain more structural information of the stable fiber, the final coarse grained configuration was backtransformed to atomistic detail. Based on this structure a 10 ns atomistic simulation was performed, which suggests that the fiber is stabilized by hydrogen bonds and water mediated hydrogen bonds. These stabilizing bonds are, however, reduced by competitive protein-water hydrogen bonds. Hence, π-stacking is suspected to play a larger role in fiber stabilization. The π-stacking of intermolecular Phe residues are found to favor a T-shaped stacking mode, while intramolecular π-stacking interactions assume a broad variety of modes from the parallel displaced mode to the T-shaped stacking mode and modes in between, with equal probability. Selected snapshots from the atomistic simulation were geometry optimized using semi-empirical quantum mechanical methods to validate the fiber stability and π-stacking configuration. An average Cα-RMSD was determined to be 2.68 Å. These findings indicate that the fiber may be used as a novel model system for the study of amyloid fibers or self-assembled conductive biowires, respectively.

Hydrolysis of Guanosine Triphosphate (GTP) by the Ras·GAP Protein Complex: Reaction Mechanism and Kinetic Scheme

Molecular mechanisms of the hydrolysis of guanosine triphosphate (GTP) to guanosine diphosphate (GDP) and inorganic phosphate (Pi) by the Ras·GAP protein complex are fully investigated by using modern modeling tools. The previously hypothesized stages of the cleavage of the phosphorus-oxygen bond in GTP and the formation of the imide form of catalytic Gln61 from Ras upon creation of Pi are confirmed by using the higher-level quantum-based calculations. The steps of the enzyme regeneration are modeled for the first time, providing a comprehensive description of the catalytic cycle. It is found that for the reaction Ras·GAP·GTP·H2O → Ras·GAP·GDP·Pi, the highest barriers correspond to the process of regeneration of the active site but not to the process of substrate cleavage. The specific shape of the energy profile is responsible for an interesting kinetic mechanism of the GTP hydrolysis. The analysis of the process using the first-passage approach and consideration of kinetic equations suggest that the overall reaction rate is a result of the balance between relatively fast transitions and low probability of states from which these transitions are taking place. Our theoretical predictions are in excellent agreement with available experimental observations on GTP hydrolysis rates.

Improving intermolecular interactions in DFTB3 using extended polarization from chemical-potential equalization

Semi-empirical quantum mechanical methods traditionally expand the electron density in a minimal, valence-only electron basis set. The minimal-basis approximation causes molecular polarization to be underestimated, and hence intermolecular interaction energies are also underestimated, especially for intermolecular interactions involving charged species. In this work, the third-order self-consistent charge density functional tight-binding method (DFTB3) is augmented with an auxiliary response density using the chemical-potential equalization (CPE) method and an empirical dispersion correction (D3). The parameters in the CPE and D3 models are fitted to high-level CCSD(T) reference interaction energies for a broad range of chemical species, as well as dipole moments calculated at the DFT level; the impact of including polarizabilities of molecules in the parameterization is also considered. Parameters for the elements H, C, N, O, and S are presented. The Root Mean Square Deviation (RMSD) interaction energy is improved from 6.07 kcal/mol to 1.49 kcal/mol for interactions with one charged species, whereas the RMSD is improved from 5.60 kcal/mol to 1.73 for a set of 9 salt bridges, compared to uncorrected DFTB3. For large water clusters and complexes that are dominated by dispersion interactions, the already satisfactory performance of the DFTB3-D3 model is retained; polarizabilities of neutral molecules are also notably improved. Overall, the CPE extension of DFTB3-D3 provides a more balanced description of different types of non-covalent interactions than Neglect of Diatomic Differential Overlap type of semi-empirical methods (e.g., PM6-D3H4) and PBE-D3 with modest basis sets.