Free Energy Methods Research Articles

Comprehensive Evaluation of End-Point Free Energy Techniques in Carboxylated-Pillar[6]arene Host-Guest Binding: IV. The QM Treatment, GB Models and the Multi-Trajectory Extension

Due to the similarity of host–guest complexes and protein–ligand and protein–protein assemblies, computational tools for protein–drug complexes are commonly applied in host–guest binding. One of the methods with the highest popularity is the end-point free energy technique, which estimates the binding affinity with gas-phase and solvation contributions extracted from simplified end-point sampling. Our series papers on a set of carboxylated-pillararene host–guest complexes have proven with solid numerical evidence that standard end-point techniques are practically useless in host–guest binding, but alterations, such as slightly increasing interior dielectric constant in post-processing calculation and shifting to the multi-trajectory realization in conformational sampling, could better the situation and pull the end-point method back to the pool of usable tools. Also, the force-field selection plays a critical role, as it determines the sampled region in the conformational space. In the current work, we continue the efforts to explore potentially promising end-point modifications in host–guest binding and further extend the sampling time to an unprecedent length. Specifically, we comprehensively benchmarked the shift from the original MM description to QM Hamiltonians in post-processing the popular single-trajectory sampling. Two critical settings in the multi-scale QM/GBSA regime are the selections of the QM Hamiltonian and the implicit-solvent model, and a scan of combinations of popular semi-empirical QM Hamiltonians and GB models is performed. The multi-scale QM/GBSA treatment is further combined with the three-trajectory sampling protocol, introducing a further advanced modification. The sampling lengths in the host–guest complex, solvated guest and solvated host ensembles are extended to 500 ns, 500 ns and 12,000 ns. As a result, the sampling quality in end-point calculations is unprecedently high, enabling us to draw conclusive pictures of investigated forms of modified end-point free energy methods. Numerical results suggest that the shift to the QM Hamiltonian does not better the situation in the popular single-trajectory regime, but noticeable improvements are observed in the three-trajectory sampling regime, especially for the DFTB/GBSA parameter combination (either DFTB2 or its third-order extension), the quality metrics of which reach an unprecedently high level and surpass existing predictions (including costly alchemical transformations) on this dataset, hinting on the applicability of the advanced three-trajectory QM/GBSA end-point modification for host–guest complexes.

Mechanical, thermodynamical, and microstructural characterization of adhesion evolution in asphalt binder-aggregate interface subjected to calcium chloride deicer

Frigid weather and freezing road conditions in cold regions result in reduced pavement surface friction, compromised vehicle maneuverability, leading to deteriorated transportation safety. Given their effectiveness and low cost, using chloride-based deicing agents such as CaCl2 has become a standard winter maintenance strategy for many highway agencies. However, coupled with freeze–thaw (F-T) cycles, CaCl2 can adversely affect the durability of asphalt pavements and cause complex distresses such as moisture-induced damage, among others. Despite this fact, the mechanisms and the extent of the damage made to the asphalt binder-aggregate interface due to use of chloride salts are not well-understood. This study was undertaken to characterize the mechanisms by which different eutectic concentrations of aqueous CaCl2 solutions and F-T cycles affect the adhesion between asphalt binder and different aggregates through a multi-scale characterization and testing program. Given the rise in use of polymer-modified asphalt binder a PG 64–34 asphalt binder and aggregates of granite, quartzite, and limestone mineralogies were evaluated. More specifically, pull-off strength of the aggregate-binder systems subjected to different aqueous concentrations of CaCl2 and F-T cycles were determined. In addition, surface free energy method as a thermodynamic-based approach was employed to determine the effect of salt concentration on the adhesion and debonding energies, and the moisture-induced damage potential of the binder-aggregate interfaces. Furthermore, atomic force microscopy was employed to evaluate the effect of salt concentration and F-T cycles on surface morphology, nanostructure, and adhesion of the asphalt binder. The applied multi-scale characterization provided invaluable information about the physio-chemical phenomena responsible for the accelerated damage in asphalt binder-aggregate interface when subjected to different aqueous concentrations of CaCl2. Gaining an accurate understanding of the damage mechanism is necessary for developing new materials and methods for combatting salt-accelerated moisture-induced damage in cold regions as a result of winter maintenance operations.

The investigation of surface free energy components and moisture sensitivity damage of asphalt mixes modified with carbon black using the sessile drop method

Asphalt pavement is vulnerable to moisture damage when water infiltrates the asphalt-aggregate interface and weakens their bond. To prevent this problem, researchers have explored the use of additives to increase resistance against moisture damage. The moisture sensitivity of asphalt mixtures can be assessed by measuring their basic properties related to the mechanism of moisture damage. One effective approach to this is using the surface free energy method, which calculates bond energies between different materials’ surfaces and identifies the most appropriate combination of asphalt and aggregate. High cohesion or adhesion bonding is crucial for the proper performance of asphalt mixtures, such as moisture damage resistance. This study aimed to investigate the effect of carbon black additive on moisture damage resistance using surface free energy methods, comparing results with the Tensile Strength Ratio (TSR) test. Modified asphalts containing 4, 8, 12, and 16% carbon black were tested, with natural base asphalt as a control sample. The study used surface free energy theory to determine moisture sensitivity parameters of asphalt mixtures by obtaining surface free energy components of the raw materials and testing the ratio of tensile strength. Results showed that adding 8% carbon black to the asphalt mixture improved moisture sensitivity, and increasing additive content up to 16% reduced it but still provided better results than samples made with pure asphalt. The study concluded that the surface free energy method can determine the suitable combination of asphalt and aggregate materials for resistance against moisture damage.

Accurately Predicting Protein pKa Values Using Nonequilibrium Alchemy.

The stability, solubility, and function of a protein depend on both its net charge and the protonation states of its individual residues. pKa is a measure of the tendency for a given residue to (de)protonate at a specific pH. Although pKa values can be resolved experimentally, theory and computation provide a compelling alternative. To this end, we assess the applicability of a nonequilibrium (NEQ) alchemical free energy method to the problem of pKa prediction. On a data set of 144 residues that span 13 proteins, we report an average unsigned error of 0.77 ± 0.09, 0.69 ± 0.09, and 0.52 ± 0.04 pK for aspartate, glutamate, and lysine, respectively. This is comparable to current state-of-the-art predictors and the accuracy recently reached using free energy perturbation methods (e.g., FEP+). Moreover, we demonstrate that our open-source, pmx-based approach can accurately resolve the pKa values of coupled residues and observe a substantial performance disparity associated with the lysine partial charges in Amber14SB/Amber99SB*-ILDN, for which an underused fix already exists.

Free energy along drug-protein binding pathways interactively sampled in virtual reality

We describe a two-step approach for combining interactive molecular dynamics in virtual reality (iMD-VR) with free energy (FE) calculation to explore the dynamics of biological processes at the molecular level. We refer to this combined approach as iMD-VR-FE. Stage one involves using a state-of-the-art ‘human-in-the-loop’ iMD-VR framework to generate a diverse range of protein–ligand unbinding pathways, benefitting from the sophistication of human spatial and chemical intuition. Stage two involves using the iMD-VR-sampled pathways as initial guesses for defining a path-based reaction coordinate from which we can obtain a corresponding free energy profile using FE methods. To investigate the performance of the method, we apply iMD-VR-FE to investigate the unbinding of a benzamidine ligand from a trypsin protein. The binding free energy calculated using iMD-VR-FE is similar for each pathway, indicating internal consistency. Moreover, the resulting free energy profiles can distinguish energetic differences between pathways corresponding to various protein–ligand conformations (e.g., helping to identify pathways that are more favourable) and enable identification of metastable states along the pathways. The two-step iMD-VR-FE approach offers an intuitive way for researchers to test hypotheses for candidate pathways in biomolecular systems, quickly obtaining both qualitative and quantitative insight.

Modern Alchemical Free Energy Methods for Drug Discovery Explained.

This Perspective provides a contextual explanation of the current state-of-the-art alchemical free energy methods and their role in drug discovery as well as highlights select emerging technologies. The narrative attempts to answer basic questions about what goes on "under the hood" in free energy simulations and provide general guidelines for how to run simulations and analyze the results. It is the hope that this work will provide a valuable introduction to students and scientists in the field.

Computing Free Energies of Fold-Switching Proteins Using MELD x MD.

Some proteins are conformational switches, able to transition between relatively different conformations. To understand what drives them requires computing the free-energy difference ΔGAB between their stable states, A and B. Molecular dynamics (MD) simulations alone are often slow because they require a reaction coordinate and must sample many transitions in between. Here, we show that modeling employing limited data (MELD) x MD on known endstates A and B is accurate and efficient because it does not require passing over barriers or knowing reaction coordinates. We validate this method on two problems: (1) it gives correct relative populations of α and β conformers for small designed chameleon sequences of protein G; and (2) it correctly predicts the conformations of the C-terminal domain (CTD) of RfaH. Free-energy methods like MELD x MD can often resolve structures that confuse machine-learning (ML) methods.

In Silico and In Vitro Potential of FDA-Approved Drugs for Antimalarial Drug Repurposing against Plasmodium Serine Hydroxymethyltransferases.

Malaria has spread in many countries, with a 12% increase in deaths after the coronavirus disease 2019 pandemic. Malaria is one of the most concerning diseases in the Greater Mekong subregion, showing increased drug-resistant rates. Serine hydroxymethyltransferase (SHMT), a key enzyme in the deoxythymidylate synthesis pathway, has been identified as a promising antimalarial drug target due to its conserved folate binding pocket. This study used a molecular docking approach to screen 2509 US Food and Drug Administration (FDA)-approved drugs against seven Plasmodium SHMT structures. Eight compounds had significantly lower binding energies than the known SHMT inhibitors pyrazolopyran(+)-86, tetrahydrofolate, and antimalarial drugs, ranging from 4 to 10 kcal/mol. Inhibition assays testing the eight compounds against Plasmodium falciparum SHMT (PfSHMT) showed that amphotericin B was a competitive inhibitor of PfSHMT with a half-maximal inhibitory concentration (IC50) of 106 ± 1 μM. Therefore, a 500 ns molecular dynamics simulation of PfSHMT/PLS/amphotericin B was performed. The backbone root-mean-square deviation of the protein-ligand complex indicated the high complex stability during simulations, supported by its radius of gyration, hydrogen-bond interactions, and number of atom contacts. The appreciable binding affinity of amphotericin B for PfSHMT was indicated by their solvated interaction energy (-11.15 ± 0.09 kcal/mol) and supported by strong ligand-protein interactions (≥80% occurrences) with its essential residues (i.e., Y78, K151, N262, F266, and V365) predicted by pharmacophore modeling and per-residue decomposition free energy methods. Therefore, our findings identify a promising new PfSHMT inhibitor, albeit with less inhibitory activity, and suggest a core structure that differs from that of previous SHMT inhibitors, thus being a rational approach for novel antimalarial drug design.

Alchemical Free Energy Workflows for the Computation of Protein-Ligand Binding Affinities.

Alchemical free energy methods can be used for the efficient computation of relative binding free energies during preclinical drug discovery stages. In recent years, this has been facilitated further by the implementation of workflows that enable non-experts to quickly and consistently set up the required simulations. Given the correct input structures, workflows handle the difficult aspects of setting up perturbations, including consistently defining the perturbable molecule, its atom mapping and topology generation, perturbation network generation, running of the simulations via different sampling methods, and analysis of the results. Different academic and commercial workflows are discussed, including FEW, FESetup, FEPrepare, CHARMM-GUI, Transformato, PMX, QLigFEP, TIES, ProFESSA, PyAutoFEP, BioSimSpace, FEP+, Flare, and Orion. These workflows differ in various aspects, such as mapping algorithms or enhanced sampling methods. Some workflows can accommodate more than one molecular dynamics (MD) engine and use external libraries for tasks. Differences between workflows can present advantages for different use cases, however a lack of interoperability of the workflows' components hinders systematic comparisons.

Design, synthesis, computational study and cytotoxic evaluation of some new quinazoline derivatives containing pyrimidine moiety

Quinazoline derivatives, as an important category of heterocyclic compounds, have received much attention for the design and development of new drugs due to their various pharmacological properties. Besides, there is a great deal of evidence showing pyrimidine analogs as anticancer agents. Thus, in the present study, for the design of new target compounds with cytotoxic activity, we focused on various quinazolinone and pyrimidine hybrids. A new series of quinazoline-pyrimidine hybrid derivatives (6a-6n) have been designed and synthesized as novel antiproliferative agents. All the synthesized compounds characterized based on their IR, NMR and Mass spectroscopic data. Antiproliferative activities of the new compounds were evaluated against three human cancer cell lines (MCF-7, A549, SW-480). The compounds were found to have appropriate potential with IC50 values ranging from 2.3 ± 5.91 to 176.5 ± 0.7 μM against the tested cell lines. Compound 6n exerted the highest antiproliferative activity with IC50 values of 5.9 ± 1.69 μM, 2.3 ± 5.91 μM and 5.65 ± 2.33 μM against A549, SW-480 and MCF-7 respectively. The results indicated that 6n could induce apoptosis in A549 cell line in a dose dependent manner and arrest in the S phase of cell cycle. Docking studies were also done to investigate the detailed binding pattern of the synthesized compounds against EGFR. Furthermore, molecular dynamic simulation and binding free energy calculation have been done to rescore initial docking pose of the synthesized compounds using ensemble-based MMGB/PBSA free energy method. According to the results, free energy calculation confirmed biological activity of compounds and also, Arg 817 and Lys 721 residues had the pivotal role in the high potency of 6n. Finally, the drug likeness and in silico ADME study were also predicted.

Exploring potential Plasmodium kinase inhibitors: a combined docking, MD and QSAR studies

Malaria is a disease caused mostly by Plasmodium falciparum, affects millions of people each year. The kinases are validated targets for malaria infection. In this study, we investigate for real and hypothetical compounds that can inhibit cyclic guanosine monophosphate (CGMP)-dependent protein kinase using molecular docking via combined similarity analysis, molecular dynamics simulations, quantitative structure activity relationship (QSAR). Using Tanimoto similarity scores, ∼8.4 million compounds were screened. Compounds that have at least 70% similarity are used in further analysis. These compounds are assessed by means of docking, MMBPSA, MMGBSA and ANI_LIE. Based on consensus of different free energy methods and docking we revealed two potential inhibitors that can be useful for treatment of malaria. Apart from screening of real compounds, we have also selected the 10 most plausible hypothetical compounds by performing QSAR. By QSAR proposed pharmacophores, we generated over 247 hypothetical compounds and among them 19 molecules with lower QSAR predicted IC50 values and high docking scores were selected for further analysis. We selected the top 10 inhibitor candidates and performed MD simulations for free energy calculations like the protocol applied for real compounds. According to the free energy calculations, we suggest 2 real (C34H29F5N8O4S and C30H27F2N7O2S2, PubChem IDs: 140564801 and 89035196, respectively) and 2 hypothetical (C23H27FN6O2S, MOL3 and C23H25FN6O2S, MOL4) compounds that can be effective inhibitors against the protein kinase of Plasmodium falciparum. Communicated by Ramaswamy H. Sarma

Thermodynamics of healing of asphalt binders with free energy

Inspired by self-healing attributes, one of the key topics in the field of pavement engineering is changing the way that asphalt materials behave for prolonged service life and environmental savings. Extensive experimental approaches, models and technologies have been developed to characterize the evolution of healing influenced by internal and external factors. The purpose of this study is to explore the local kinetics of asphalt binders that energetically favor autonomous self-healing as well as its mechanism using free energy method. With consideration of driving forces including internal stress and interfacial forces of attraction (or surface energy) responsible for energy redistribution during the rest period, the healing of an asphalt binder was first modeled as a function of thermodynamic parameters. New healing indexes were then defined by the total free energy change for healing (ΔGH) and fast-rate response ratio (RH) of initial healing in duration to the total rest time. Thirdly, a fatigue-healing test was performed on #70 asphalt binder and SBS-modified asphalt binder to examine the potential and reactivity of healing affected by different damage levels. Finally, the effectiveness of proposed thermodynamic healing indexes in describing healing ability was discussed by considering their relations with properties measured by mechanical analysis. Results show that more negative ΔGH indicates a greater tendency to spontaneous healing while less negative ΔGH suggests a poor healing potential as crack propagates to a certain length. In contrast to #70 asphalt binder, SBS-modified asphalt binder exhibits a superior healing potential on a shorter crack length scale. RH gradually increases with the increase of crack length for both #70 asphalt binder and SBS-modified asphalt binder, implying limited mobility of molecules. By contrast with mechanistic healing indexes assessing properties before and after healing (i.e. the degree of healing), the thermodynamic healing indexes provide a new vision regarding healing potential, which reveals the true ability of a material to autonomously repair damaged areas.

Identification of oil palm cis-regulatory elements based on DNA free energy and single nucleotide polymorphism density

Transcription control through cis-regulatory elements (CREs) is one of important regulators of gene expression. This study aimed to identify the location of CREs in oil palm (Elaeis guineensis Jacq.) using the combination of DNA free energy and single nucleotide polymorphism (SNP) density approaches. Promoter region sequences were extracted oil palm genome spanning from 1500 nucleotides (nt) upstream to 1000 nt downstream of every annotated transcription start sites (TSS). Free energy profiles of each promoter region were calculated using PromPredict software. Raw reads from the deep sequencing of 59 oil palm origins were used to calculate SNP density of each promoter region. The result showed that the average free energy (AFE) on the upstream region of TSS is about 1.5 kcal/mol higher compared to the downstream region. Using DNA free energy method, 16,281 regions of CREs were predicted. Most of predicted CREs was located between 1 and 500 nt upstream of TSS. Anti-correlation pattern between free energy and SNP density was observed on the predicted regions of CREs. This anti-correlated pattern was also observed on an experimentally determined promoter of the oil palm metallothionein gene, EgMSP1. Considering the increasing use of promoter information on plant biotechnology, an easy and accurate promoter prediction using the combination of free energy and SNP density method could be recommended.

Uni-GBSA: an open-source and web-based automatic workflow to perform MM/GB(PB)SA calculations for virtual screening.

Binding free energy calculation of a ligand to a protein receptor is a fundamental objective in drug discovery. Molecular mechanics/Generalized-Born (Poisson-Boltzmann) surface area (MM/GB(PB)SA) is one of the most popular methods for binding free energy calculations. It is more accurate than most scoring functions and more computationally efficient than alchemical free energy methods. Several open-source tools for performing MM/GB(PB)SA calculations have been developed, but they have limitations and high entry barriers to users. Here, we introduce Uni-GBSA, a user-friendly automatic workflow to perform MM/GB(PB)SA calculations, which can perform topology preparation, structure optimization, binding free energy calculation and parameter scanning for MM/GB(PB)SA calculations. It also offers a batch mode that evaluates thousands of molecules against one protein target in parallel for efficient application in virtual screening. The default parameters are selected after systematic testing on the PDBBind-2011 refined dataset. In our case studies, Uni-GBSA produced a satisfactory correlation with the experimental binding affinities and outperformed AutoDock Vina in molecular enrichment. Uni-GBSA is available as an open-source package at https://github.com/dptech-corp/Uni-GBSA. It can also be accessed for virtual screening from the Hermite web platform at https://hermite.dp.tech. A free Uni-GBSA web server of a lab version is available at https://labs.dp.tech/projects/uni-gbsa/. This increases user-friendliness because the web server frees users from package installations and provides users with validated workflows for input data and parameter settings, cloud computing resources for efficient job completions, a user-friendly interface and professional support and maintenance.

Bardeen black hole thermodynamics from topological perspective

In this paper, we use the generalized off-shell Helmholtz free energy method to explore the thermodynamic properties of Bardeen black holes (BD BHs) from a topological perspective based on Duan’s topological current ϕ-mapping. We consider various structures of BD BHs, including regular BD-AdS BHs, BD-AdS BHs in Kiselev’s model of quintessence, BD BHs in massive gravity (MG), and BD BHs in 4D Einstein–Gauss–Bonnet (EGB) gravity. We demonstrate that these BHs have one topological classification (TC), i.e., TC is +1 for all BHs considered, and the addition of MG or GB terms, etc., does not change the topological numbers.

PyCHARMM: Embedding CHARMM Functionality in a Python Framework.

CHARMM is rich in methodology and functionality as one of the first programs addressing problems of molecular dynamics and modeling of biological macromolecules and their partners, e.g., small molecule ligands. When combined with the highly developed CHARMM parameters for proteins, nucleic acids, small molecules, lipids, sugars, and other biologically relevant building blocks, and the versatile CHARMM scripting language, CHARMM has been a trendsetting platform for modeling studies of biological macromolecules. To further enhance the utility of accessing and using CHARMM functionality in increasingly complex workflows associated with modeling biological systems, we introduce pyCHARMM, Python bindings, functions, and modules to complement and extend the extensive set of modeling tools and methods already available in CHARMM. These include access to CHARMM function-generated variables associated with the system (psf), coordinates, velocities and forces, atom selection variables, and force field related parameters. The ability to augment CHARMM forces and energies with energy terms or methods derived from machine learning or other sources, written in Python, CUDA, or OpenCL and expressed as Python callable routines is introduced together with analogous functions callable during dynamics calculations. Integration of Python-based graphical engines for visualization of simulation models and results is also accessible. Loosely coupled parallelism is available for workflows such as free energy calculations, using MBAR/TI approaches or high-throughput multisite λ-dynamics (MSλD) free energy methods, string path optimization calculations, replica exchange, and molecular docking with a new Python-based CDOCKER module. CHARMM accelerated platform kernels through the CHARMM/OpenMM API, CHARMM/DOMDEC, and CHARMM/BLaDE API are also readily integrated into this Python framework. We anticipate that pyCHARMM will be a robust platform for the development of comprehensive and complex workflows utilizing Python and its extensive functionality as well as an optimal platform for users to learn molecular modeling methods and practices within a Python-friendly environment such as Jupyter Notebooks.

Understanding the sources of error in MBAR through asymptotic analysis.

Many sampling strategies commonly used in molecular dynamics, such as umbrella sampling and alchemical free energy methods, involve sampling from multiple states. The Multistate Bennett Acceptance Ratio (MBAR) formalism is a widely used way of recombining the resulting data. However, the error of the MBAR estimator is not well-understood: previous error analyses of MBAR assumed independent samples. In this work, we derive a central limit theorem for MBAR estimates in the presence of correlated data, further justifying the use of MBAR in practical applications. Moreover, our central limit theorem yields an estimate of the error that can be decomposed into contributions from the individual Markov chains used to sample the states. This gives additional insight into how sampling in each state affects the overall error. We demonstrate our error estimator on an umbrella sampling calculation of the free energy of isomerization of the alanine dipeptide and an alchemical calculation of the hydration free energy of methane. Our numerical results demonstrate that the time required for the Markov chain to decorrelate in individual states can contribute considerably to the total MBAR error, highlighting the importance of accurately addressing the effect of sample correlation.

Long Time Scale Ensemble Methods in Molecular Dynamics: Ligand-Protein Interactions and Allostery in SARS-CoV-2 Targets.

We subject a series of five protein-ligand systems which contain important SARS-CoV-2 targets, 3-chymotrypsin-like protease (3CLPro), papain-like protease, and adenosine ribose phosphatase, to long time scale and adaptive sampling molecular dynamics simulations. By performing ensembles of ten or twelve 10 μs simulations for each system, we accurately and reproducibly determine ligand binding sites, both crystallographically resolved and otherwise, thereby discovering binding sites that can be exploited for drug discovery. We also report robust, ensemble-based observation of conformational changes that occur at the main binding site of 3CLPro due to the presence of another ligand at an allosteric binding site explaining the underlying cascade of events responsible for its inhibitory effect. Using our simulations, we have discovered a novel allosteric mechanism of inhibition for a ligand known to bind only at the substrate binding site. Due to the chaotic nature of molecular dynamics trajectories, regardless of their temporal duration individual trajectories do not allow for accurate or reproducible elucidation of macroscopic expectation values. Unprecedentedly at this time scale, we compare the statistical distribution of protein-ligand contact frequencies for these ten/twelve 10 μs trajectories and find that over 90% of trajectories have significantly different contact frequency distributions. Furthermore, using a direct binding free energy calculation protocol, we determine the ligand binding free energies for each of the identified sites using long time scale simulations. The free energies differ by 0.77 to 7.26 kcal/mol across individual trajectories depending on the binding site and the system. We show that, although this is the standard way such quantities are currently reported at long time scale, individual simulations do not yield reliable free energies. Ensembles of independent trajectories are necessary to overcome the aleatoric uncertainty in order to obtain statistically meaningful and reproducible results. Finally, we compare the application of different free energy methods to these systems and discuss their advantages and disadvantages. Our findings here are generally applicable to all molecular dynamics based applications and not confined to the free energy methods used in this study.

Best Practices on QM/MM Simulations of Biological Systems.

During the second half of the 20th century, following structural biology hallmark works on DNA and proteins, biochemists shifted their questions from "what does this molecule look like?" to "how does this process work?". Prompted by the theoretical and practical developments in computational chemistry, this led to the emergence of biomolecular simulations and, along with the 2013 Nobel Prize in Chemistry, to the development of hybrid QM/MM methods. QM/MM methods are necessary whenever the problem we want to address involves chemical reactivity and/or a change in the system's electronic structure, with archetypal examples being the studies of an enzyme's reaction mechanism and a metalloprotein's active site. In the last decades QM/MM methods have seen an increasing adoption driven by their incorporation in widely used biomolecular simulation software. However, properly setting up a QM/MM simulation is not an easy task, and several issues need to be properly addressed to obtain meaningful results. In the present work, we describe both the theoretical concepts and practical issues that need to be considered when performing QM/MM simulations. We start with a brief historical perspective on the development of these methods and describe when and why QM/MM methods are mandatory. Then we show how to properly select and analyze the performance of the QM level of theory, the QM system size, and the position and type of the boundaries. We show the relevance of performing prior QM model system (or QM cluster) calculations in a vacuum and how to use the corresponding results to adequately calibrate those derived from QM/MM. We also discuss how to prepare the starting structure and how to select an adequate simulation strategy, including those based on geometry optimizations as well as free energy methods. In particular, we focus on the determination of free energy profiles using multiple steered molecular dynamics (MSMD) combined with Jarzynski's equation. Finally, we describe the results for two illustrative and complementary examples: the reaction performed by chorismate mutase and the study of ligand binding to hemoglobins. Overall, we provide many practical recommendations (or shortcuts) together with important conceptualizations that we hope will encourage more and more researchers to incorporate QM/MM studies into their research projects.

MELD-Bracket Ranks Binding Affinities of Diverse Sets of Ligands.

Affinity ranking of structurally diverse small-molecule ligands is a challenging problem with important applications in structure-based drug discovery. Absolute binding free energy methods can model diverse ligands, but the high computational cost of the current methods limits application to data sets with few ligands. We recently developed MELD-Bracket, a Molecular Dynamics method for efficient affinity ranking of ligands [ JCTC 2022, 18 (1), 374-379]. It utilizes a Bayesian framework to guide sampling to relevant regions of phase space, and it couples this with a bracket-like competition on a pool of ligands. Here we find that 6-competitor MELD-Bracket can rank dozens of diverse ligands that have low structural similarity and different net charges. We benchmark it on four protein systems─PTB1B, Tyk2, BACE, and JAK3─having varied modes of interactions. We also validated 8-competitor and 12-competitor protocols. The MELD-Bracket protocols presented here may have the appropriate balance of accuracy and computational efficiency to be suitable for ranking diverse ligands from typical drug discovery campaigns.