QM/MM Free Energy Calculations of Long-Range Biological Protonation Dynamics by Adaptive and Focused Sampling.

Water-mediated proton transfer reactions are central for catalytic processes in a wide range of biochemical systems, ranging from biological energy conversion to chemical transformations in the metabolism. Yet, the accurate computational treatment of such complex biochemical reactions is highly challenging and requires the application of multiscale methods, in particular hybrid quantum/classical (QM/MM) approaches combined with free energy simulations. Here, we combine the unique exploration power of new advanced sampling methods with density functional theory (DFT)-based QM/MM free energy methods for multiscale simulations of long-range protonation dynamics in biological systems. In this regard, we show that combining multiple walkers/well-tempered metadynamics with an extended system adaptive biasing force method (MWE) provides a powerful approach for exploration of water-mediated proton transfer reactions in complex biochemical systems. We compare and combine the MWE method also with QM/MM umbrella sampling and explore the sampling of the free energy landscape with both geometric (linear combination of proton transfer distances) and physical (center of excess charge) reaction coordinates and show how these affect the convergence of the potential of mean force (PMF) and the activation free energy. We find that the QM/MM-MWE method can efficiently explore both direct and water-mediated proton transfer pathways together with forward and reverse hole transfer mechanisms in the highly complex proton channel of respiratory Complex I, while the QM/MM-US approach shows a systematic convergence of selected long-range proton transfer pathways. In this regard, we show that the PMF along multiple proton transfer pathways is recovered by combining the strengths of both approaches in a QM/MM-MWE/focused US (FUS) scheme and reveals new mechanistic insight into the proton transfer principles of Complex I. Our findings provide a promising basis for the quantitative multiscale simulations of long-range proton transfer reactions in biological systems.

Systematic QM/MM Study for Predicting 31P NMR Chemical Shifts of Adenosine Nucleotides in Solution and Stages of ATP Hydrolysis in a Protein Environment.

NMR (nuclear magnetic resonance) spectroscopy allows for important atomistic insights into the structure and dynamics of biological macromolecules; however, reliable assignments of experimental spectra are often difficult. Herein, quantum mechanical/molecular mechanical (QM/MM) calculations can provide crucial support. A major problem for the simulations is that experimental NMR signals are time-averaged over much longer time scales, and since computed chemical shifts are highly sensitive to local changes in the electronic and structural environment, sufficiently large averages over representative structural ensembles are essential. This entails high computational demands for reliable simulations. For NMR measurements in biological systems, a nucleus of major interest is 31P since it is both highly present (e.g., in nucleic acids) and easily observable. The focus of our present study is to develop a robust and computationally cost-efficient framework for simulating 31P NMR chemical shifts of nucleotides. We apply this scheme to study the different stages of the ATP hydrolysis reaction catalyzed by p97. Our methodology is based on MM molecular dynamics (MM-MD) sampling, followed by QM/MM structure optimizations and NMR calculations. Overall, our study is one of the most comprehensive QM-based 31P studies in a protein environment and the first to provide computed NMR chemical shifts for multiple nucleotide states in a protein environment. This study sheds light on a process that is challenging to probe experimentally and aims to bridge the gap between measured and calculated NMR spectroscopic properties.

Entropy and Energy Profiles of Chemical Reactions.

The description of chemical processes at the molecular level is often facilitated by the use of reaction coordinates or collective variables (CVs). The CV measures the progress of the reaction and allows the construction of profiles that track how specific properties evolve as the reaction progresses. Whereas CVs are routinely used, especially alongside enhanced sampling techniques, the links among reaction profiles, thermodynamic state functions, and reaction rate constants are not rigorously exploited. Here, we report a unified treatment of such reaction profiles. Tractable expressions are derived for the free-energy, internal-energy, and entropy profiles as functions of only the CV. We demonstrate the ability of this treatment to extract quantitative insight from the entropy and internal-energy profiles of various real-world physicochemical processes, including intramolecular organic reactions, ionic transport in superionic electrolytes, and molecular transport in nanoporous materials.

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.

From free-energy profiles to activation free energies.

Given a chemical reaction going from reactant (R) to the product (P) on a potential energy surface (PES) and a collective variable (CV) discriminating between R and P, we define the free-energy profile (FEP) as the logarithm of the marginal Boltzmann distribution of the CV. This FEP is not a true free energy. Nevertheless, it is common to treat the FEP as the "free-energy" analog of the minimum potential energy path and to take the activation free energy, ΔFRP ‡, as the difference between the maximum at the transition state and the minimum at R. We show that this approximation can result in large errors. The FEP depends on the CV and is, therefore, not unique. For the same reaction, different discriminating CVs can yield different ΔFRP ‡. We derive an exact expression for the activation free energy that avoids this ambiguity. We find ΔFRP ‡ to be a combination of the probability of the system being in the reactant state, the probability density on the dividing surface, and the thermal de Broglie wavelength associated with the transition. We apply our formalism to simple analytic models and realistic chemical systems and show that the FEP-based approximation applies only at low temperatures for CVs with a small effective mass. Most chemical reactions occur on complex, high-dimensional PES that cannot be treated analytically and pose the added challenge of choosing a good CV. We study the influence of that choice and find that, while the reaction free energy is largely unaffected, ΔFRP ‡ is quite sensitive.

Accelerating Hybrid Density Functional Theory Molecular Dynamics Simulations by Seminumerical Integration, Resolution-of-the-Identity Approximation, and Graphics Processing Units.

The computationally very demanding evaluation of the 4-center-2-electron (4c2e) integrals and their respective integral derivatives typically represents the major bottleneck within hybrid Kohn-Sham density functional theory molecular dynamics simulations. Building upon our previous works on seminumerical exact-exchange (sn-LinK) [Laqua, H., Thompsons, T. H., Kussmann, J., Ochsenfeld, C., J. Chem. Theory Comput. 2020, 16, 1465] and resolution-of-the-identity Coulomb (RI-J) [Kussmann, J., Laqua, H., Ochsenfeld, C., J. Chem. Theory Comput. 2021, 17, 1512], the expensive 4c2e integral evaluation can be avoided entirely, resulting in a highly efficient electronic structure theory method, allowing for fast ab initio molecular dynamics (AIMD) simulations even with large basis sets. Moreover, we propose to combine the final self-consistent field (SCF) step with the subsequent nuclear forces evaluation, providing the forces at virtually no additional cost after a converged SCF calculation, reducing the total runtime of an AIMD simulation by about another 25%. In addition, multiple independent MD trajectories can be computed concurrently on a single node, leading to a greatly increased utilization of the available hardware─especially when combined with graphics processing unit acceleration─improving the overall throughput by up to another 5 times in this way. With all of those optimizations combined, our proposed method provides nearly 3 orders of magnitude faster execution times than traditional 4c2e integral-based methods. To demonstrate the practical utility of the approach, quantum-mechanical/molecular-mechanical dynamics simulations on double-stranded DNA were performed, investigating the relative hydrogen bond strength between adenine-thymine and guanine-cytosine base pairs. In addition, this illustrative application also contains a general accuracy assessment of the introduced approximations (integration grids, resolution-of-the-identity) within AIMD simulations, serving as a protocol on how to apply these new methods to practical problems.

Development of hetero-triaryls as a new chemotype for subtype-selective and potent Sirt5 inhibition.

In contrast to other sirtuins (NAD+-dependent class III lysine deacylases), inhibition of Sirt5 is poorly investigated, yet. Our present work is based on the recently identified Sirt5 inhibitor balsalazide, an approved drug with negligible bioavailability after oral administration. After gaining first insights into its structure-activity relationship in previous work, we were able to now develop heteroaryl-triaryls as a novel chemotype of drug-like, potent and subtype-selective Sirt5 inhibitors. The unfavourable azo group of the lead structure was modified in a systematic and comprehensive manner, leading us to a few open-chained and, most importantly, five-membered heteroaromatic substitutes (isoxazole CG_209, triazole CG_220, pyrazole CG_232) with very encouraging in vitro activities (IC50 on Sirt5 in the low micromolar range, <10 µM). These advanced inhibitors were free of cytotoxicity and showed favourable pharmacokinetic properties, as confirmed by permeability into mitochondria using live cell imaging experiments. Furthermore, results from calculations of the relative free binding affinities of the analogues compared to balsalazide as reference compound agreed well with the trends for inhibitory activities obtained in the in vitro experiments. Therefore, this method can be used to predict the affinity of closely related future potential Sirt5 inhibitors. Encouraged by our findings, we employed chemoproteomic selectivity profiling to confirm Sirt5 as main target of balsalazide and one of its improved analogues. An immobilised balsalazide-analogue specifically pulled down Sirt5 from whole cell lysates and competition experiments identified glutaryl-CoA dehydrogenase (GCDH) and nucleotide diphosphate kinase (NME4) as potential off-targets, once again confirming the selectivity of the novel balsalazide-derived Sirt5 inhibitors. In summary, a combination of targeted chemical synthesis, biological work, and computational studies led to a new generation of tailored Sirt5 inhibitors, which represent valuable chemical tools for the investigation of the physiological role of Sirt5, but could also serve as advanced lead structures for drug candidates for systemic use.

Statistically optimal analysis of the extended-system adaptive biasing force (eABF) method.

The extended-system adaptive biasing force (eABF) method and its newer variants offer rapid exploration of the configuration space of chemical systems. Instead of directly applying the ABF bias to collective variables, they are harmonically coupled to fictitious particles, which separates the problem of enhanced sampling from that of free energy estimation. The prevalent analysis method to obtain the potential of mean force (PMF) from eABF is thermodynamic integration. However, besides the PMF, most information is lost as the unbiased probability of visited configurations is never recovered. In this contribution, we show how statistical weights of individual frames can be computed using the Multistate Bennett's Acceptance Ratio (MBAR), putting the post-processing of eABF on one level with other frequently used sampling methods. In addition, we apply this formalism to the prediction of nuclear magnetic resonance shieldings, which are very sensitive to molecular geometries and often require extensive sampling. The results show that the combination of enhanced sampling by means of extended-system dynamics with the MBAR estimator is a highly useful tool for the calculation of ensemble properties. Furthermore, the extension of the presented scheme to the recently published Gaussian-accelerated molecular dynamics eABF hybrid is straightforward and approximation free.

Exponential averaging versus umbrella sampling for computing the QM/MM free energy barrier of the initial step of the desuccinylation reaction catalyzed by sirtuin 5.

The computational characterization of enzymatic reactions poses a great challenge which arises from the high dimensional and often rough potential energy surfaces commonly explored by static QM/MM methods such as adiabatic mapping (AM). The present study highlights the difficulties in estimating free energy barriers via exponential averaging over AM pathways. Based on our previous study [von der Esch et al., J. Chem. Theory Comput., 2019, 15, 6660-6667], where we analyzed the first reaction step of the desuccinylation reaction catalyzed by human sirtuin 5 (SIRT5) by means of QM/MM adiabatic mapping and machine learning, we use, here, umbrella sampling to compute the free energy profile of the initial reaction step. The computational investigations show that the initial step of the desuccinylation reaction proceeds via an SN2-type reaction mechanism in SIRT5, suggesting that the first step of the deacylation reactions catalyzed by sirtuins is highly conserved. In addition, the direct comparison of the extrapolated free energy barrier from minimal energy paths and the computed free energy path from umbrella sampling further underlines the importance of extensive sampling.

How to obtain reaction free energies from free-energy profiles.

For chemical reactions that occur via the rearrangement of atoms from a configuration about one minimum (reactant, R) of the potential energy surface (PES) to a configuration about another minimum (product, P), an exact relation between the Helmholtz reaction free energy (ΔFRP) and the free-energy profile (FEP) can be derived. Since the FEP assumes a form similar to that of the PES along the minimum energy path between R and P, there is an unfortunate tendency to regard the FEP as the "free-energy" analog of the minimum energy path and consequently to equate ΔFRP to the difference between the values of the FEP at the minima corresponding to R and P. Analytic treatments of one- and two-dimensional models are presented that show how this mistaken idea leads to errors. In effect, treating the FEP by analogy with the minimum energy path neglects the role of entropy. The FEP is a function of a collective variable (CV), which must be chosen to describe the course of the rearrangement consistently with the exact relation between ΔFRP and the FEP. For large systems of common interest, the PES is often so complex that a straightforward way of choosing a CV is lacking. Consequently, one is forced to make an educated guess. A criterion for judging the quality of the guess is proposed and applied to a two-dimensional model.

Identification of the subtype-selective Sirt5 inhibitor balsalazide through systematic SAR analysis and rationalization via theoretical investigations.

We report here an extensive structure-activity relationship study of balsalazide, which was previously identified in a high-throughput screening as an inhibitor of Sirt5. To get a closer understanding why this compound is able to inhibit Sirt5, we initially performed docking experiments comparing the binding mode of a succinylated peptide as the natural substrate and balsalazide with Sirt5 in the presence of NAD+. Based on the evidence gathered here, we designed and synthesized 13 analogues of balsalazide, in which single functional groups were either deleted or slightly altered to investigate which of them are mandatory for high inhibitory activity. Our study confirms that balsalazide with all its given functional groups is an inhibitor of Sirt5 in the low micromolar concentration range and structural modifications presented in this study did not increase potency. While changes on the N-aroyl-ß-alanine side chain eliminated potency, the introduction of a truncated salicylic acid part minimally altered potency. Calculations of the associated reaction paths showed that the inhibition potency is very likely dominated by the stability of the inhibitor-enzyme complex and not the type of inhibition (covalent vs. non-covalent). Further in-vitro characterization in a trypsin coupled assay determined that the tested inhibitors showed no competition towards NAD+ or the synthetic substrate analogue ZKsA. In addition, investigations for subtype selectivity revealed that balsalazide is a subtype-selective Sirt5 inhibitor, and our initial SAR and docking studies pave the way for further optimization.

Predicting 19 F†NMR Chemical Shifts: A Combined Computational and Experimental Study of a Trypanosomal Oxidoreductase-Inhibitor Complex.

The absence of fluorine from most biomolecules renders it an excellent probe for NMR spectroscopy to monitor inhibitor-protein interactions. However, predicting the binding mode of a fluorinated ligand from a chemical shift (or vice versa) has been challenging due to the high electron density of the fluorine atom. Nonetheless, reliable 19 F chemical-shift predictions to deduce ligand-binding modes hold great potential for in silico drug design. Herein, we present a systematic QM/MM study to predict the 19 F NMR chemical shifts of a covalently bound fluorinated inhibitor to the essential oxidoreductase tryparedoxin (Tpx) from African trypanosomes, the causative agent of African sleeping sickness. We include many protein-inhibitor conformations as well as monomeric and dimeric inhibitor-protein complexes, thus rendering it the largest computational study on chemical shifts of 19 F nuclei in a biological context to date. Our predicted shifts agree well with those obtained experimentally and pave the way for future work in this area.

Why Proline? Influence of Ring-Size on the Collagen Triple Helix.

The effect of four- and six-membered ring-size analogs (azetidine- and piperidine-2-carboxylic acid, H-Aze-OH and H-Pip-OH) of proline on the stability of the collagen triple helix was examined. Computational and nuclear magnetic resonance spectroscopic studies with model compounds and thermal denaturation experiments with collagen peptides showed that the ring-size analogs destabilize the triple helix to a similar extent by either mismatching backbone dihedral angles ϕ and ψ (Pip) or by an unfavorable trans/cis amide bond ratio (Aze).

Finding Reactive Configurations: A Machine Learning Approach for Estimating Energy Barriers Applied to Sirtuin 5.

Sirtuin 5 is a class III histone deacetylase that, unlike its classification, mainly catalyzes desuccinylation and demanoylation reactions. It is an interesting drug target that we use here to test new ideas for calculating reaction pathways of large molecular systems such as enzymes. A major issue with most schemes (e.g., adiabatic mapping) is that the resulting activation barrier height heavily depends on the chosen educt conformation. This makes the selection of the initial structure decisive for the success of the characterization. Here, we apply machine learning to a large number of molecular dynamics frames and potential energy barriers obtained by quantum mechanics/molecular mechanics calculations in order to identify (1) suitable start-conformations for reaction path calculations and (2) structural features relevant for the first step of the desuccinylation reaction catalyzed by Sirtuin 5. The latter generally aids the understanding of reaction mechanisms and important interactions in active centers. Using our novel approach, we found eleven key features that govern the reactivity. We were able to estimate reaction barriers with a mean absolute error of 3.6 kcal/mol and identified reactive configurations.

Reaction Mechanism for the N-Glycosidic Bond Cleavage of 5-Formylcytosine by Thymine DNA Glycosylase.

Thymine DNA glycosylase (TDG) initiates the base excision repair mechanism for the deamination and oxidation products of cytosine and 5-methylcytosine. This enzyme has a key role in epigenetic regulation, and its catalytic inactivation results in, e.g., mice embryo lethality. Here, we employ molecular dynamics simulations and quantum mechanics/molecular mechanics calculations to investigate the reaction mechanism of the TDG-catalyzed N-glycosidic bond hydrolysis of the modified base 5-formylcytosine. Our results reveal a reaction pathway, which in its first step features a reorganization of the substrate that lowers the barrier height for the subsequent C1'-N1 bond dissociation. The suggested mechanism is consistent with the experimental data, as it is not acid-catalyzed and proceeds through an oxocarbenium-like transition state. It also provides insights into the catalytic roles of the Thr197 and Asn140 residues.

Calculating free energies from the vibrational density of states function: Validation and critical assessment.

We explore and show the usefulness of the density of states function for computing vibrational free energies and free energy differences between small systems. Therefore, we compare this density of states integration method (DSI) to more established schemes such as Bennett's Acceptance Ratio method (BAR), the Normal Mode Analysis (NMA), and the Quasiharmonic Analysis (QHA). The strengths and shortcomings of all methods are highlighted with three numerical examples. Furthermore, the free energy of the ionization of ammonia and the mutation from serine to cysteine are computed using extensive ab initio molecular dynamics simulations. We conclude that DSI improves upon the other frequency-based methods (NMA and QHA) regarding the treatment of anharmonicity and yielding results comparable to BAR in all cases without the need for alchemical transformations. Low-frequency modes lead to larger errors indicating that long simulation times might be required for larger systems. In addition, we introduce the use of DSI for the localization of the vibrational free energy to specific atoms or residues, leading to insights into the underlying process, a unique feature that is only offered by this method.

Identifying Free Energy Hot-Spots in Molecular Transformations.

The free energy is one of the central quantities in material and natural sciences. While being well-established, e.g., in drug design or catalyst optimization, computational methods lack a straightforward way to gain deeper insights into the calculated free energy, and thus the underlying chemical or physical processes. Here, we present a generally applicable, spectrum-based ansatz that tackles this shortcoming by identifying contributions from specific atoms or groups to the vibrational free energy. We illustrate this in studies of the bromodomain-inhibitor binding and the anomeric effect in glucose providing quantitative evidence in line with chemical intuition in both cases. For the latter example we also report an experimental infrared spectrum and find excellent agreement with our simulated spectra.

Chemically Realistic Tetrahedral Lattice Models for Polymer Chains: Application to Polyethylene Oxide.

To speed up the generation of an ensemble of poly(ethylene oxide) (PEO) polymer chains in solution, a tetrahedral lattice model possessing the appropriate bond angles is used. The distance between noncovalently bonded atoms is maintained at realistic values by generating chains with an enhanced degree of self-avoidance by a very efficient Monte Carlo (MC) algorithm. Potential energy parameters characterizing this lattice model are adjusted so as to mimic realistic PEO polymer chains in water simulated by molecular dynamics (MD), which serves as a benchmark. The MD data show that PEO chains have a fractal dimension of about two, in contrast to self-avoiding walk lattice models, which exhibit the fractal dimension of 1.7. The potential energy accounts for a mild hydrophobic effect (HYEF) of PEO and for a proper setting of the distribution between trans and gauche conformers. The potential energy parameters are determined by matching the Flory radius, the radius of gyration, and the fraction of trans torsion angles in the chain. A gratifying result is the excellent agreement of the pair distribution function and the angular correlation for the lattice model with the benchmark distribution. The lattice model allows for the precise computation of the torsional entropy of the chain. The generation of polymer conformations of the adjusted lattice model is at least 2 orders of magnitude more efficient than MD simulations of the PEO chain in explicit water. This method of generating chain conformations on a tetrahedral lattice can also be applied to other types of polymers with appropriate adjustment of the potential energy function. The efficient MC algorithm for generating chain conformations on a tetrahedral lattice is available for download at https://github.com/Roulattice/Roulattice .

Low dehydroepiandrosterone and ischemic heart disease in middle-aged men: prospective results from the Massachusetts Male Aging Study.

The adrenal steroid dehydroepiandrosterone (DHEA) and its sulfate (DHEAS) have been characterized as "protective" against ischemic heart disease (IHD), especially in men, on the basis of sparse epidemiologic evidence. The authors used data from the Massachusetts Male Aging Study, a random sample prospective study of 1,709 men aged 40-70 years at baseline, to test whether serum levels of DHEA or DHEAS could predict incident IHD over a 9-year interval. At baseline (1987-1989) and follow-up (1995-1997), an interviewer-phlebotomist visited each subject in his home to obtain comprehensive health information, body measurements, and blood samples for hormone and lipid analysis. Incident IHD between baseline and follow-up was ascertained from hospital records and death registries, supplemented by self-report and evidence of medication. In the analysis sample of 1,167 men, those with serum DHEAS in the lowest quartile at baseline (<1.6 microg/ml) were significantly more likely to incur IHD by follow-up (adjusted odds ratio = 1.60, 95 percent confidence interval: 1.07, 2.39; p = 0.02), independently of a comprehensive set of known risk factors including age, obesity, diabetes, hypertension, smoking, serum lipids, alcohol intake, and physical activity. Low serum DHEA was similarly predictive. These results confirm prior evidence that low DHEA and DHEAS can predict IHD in men.

Measurement of erectile dysfunction in population-based studies: the use of a single question self-assessment in the Massachusetts Male Aging Study.

A concise, reliable means of assessing erectile dysfunction (ED) in large, multidisciplinary population-based studies is needed. A single, direct question for self-assessed ED was assessed in the population-based sample of the Massachusetts Male Aging Study (MMAS). Of the 1156 respondents to the 1995-97 MMAS follow-up evaluation, 505 were randomly selected to complete either the International Index of Erectile Function (IIEF) (n = 254), or the Brief Male Sexual Function Inventory (BMSFI) (n = 251), in addition to the single question self-assessment. The proportion not classified due to missing data was MMAS-9%, BMSFI-8%, and IIEF-18%. The single question correlated well with these other measures (r = 0.71-0.78, P < 0.001). Prevalence was similar to that based on the IIEF, agreement was moderate (kappa = 0.56-0.58), and associations with previously identified risk factors were similar for each classification. Thus, the MMAS single question may be a practical tool for population-based studies where detailed clinical measures of ED are impractical.