Lambda-ABF: Simplified, Portable, Accurate, and Cost-Effective Alchemical Free-Energy Computation.

We introduce the lambda-Adaptive Biasing Force (lambda-ABF) method for the computation of alchemical free-energy differences. We propose a software implementation and showcase it on biomolecular systems. The method arises from coupling multiple-walker adaptive biasing force with λ-dynamics. The sampling of the alchemical variable is continuous and converges toward a uniform distribution, making manual optimization of the λ schedule unnecessary. Contrary to most other approaches, alchemical free-energy estimates are obtained immediately without any postprocessing. Free diffusion of λ improves orthogonal relaxation compared to fixed-λ thermodynamic integration or free-energy perturbation. Furthermore, multiple walkers provide generic orthogonal space coverage with minimal user input and negligible computational overhead. We show that our high-performance implementations coupling the Colvars library with NAMD and Tinker-HP can address real-world cases including ligand-receptor binding with both fixed-charge and polarizable models, with a demonstrably richer sampling than fixed-λ methods. The implementation is fully open-source, publicly available, and readily usable by practitioners of current alchemical methods. Thanks to the portable Colvars library, lambda-ABF presents a unified user interface regardless of the back-end (NAMD, Tinker-HP, or any software to be interfaced in the future), sparing users the effort of learning multiple interfaces. Finally, the Colvars Dashboard extension of the visual molecular dynamics (VMD) software provides an interactive monitoring and diagnostic tool for lambda-ABF simulations.

Alchemical Free Energy Calculations on Membrane-Associated Proteins.

Membrane proteins have diverse functions within cells and are well-established drug targets. The advances in membrane protein structural biology have revealed drug and lipid binding sites on membrane proteins, while computational methods such as molecular simulations can resolve the thermodynamic basis of these interactions. Particularly, alchemical free energy calculations have shown promise in the calculation of reliable and reproducible binding free energies of protein-ligand and protein-lipid complexes in membrane-associated systems. In this review, we present an overview of representative alchemical free energy studies on G-protein-coupled receptors, ion channels, transporters as well as protein-lipid interactions, with emphasis on best practices and critical aspects of running these simulations. Additionally, we analyze challenges and successes when running alchemical free energy calculations on membrane-associated proteins. Finally, we highlight the value of alchemical free energy calculations calculations in drug discovery and their applicability in the pharmaceutical industry.

Molecular determinants of inhibition of UCP1-mediated respiratory uncoupling.

Brown adipose tissue expresses uncoupling protein 1 (UCP1), which dissipates energy as heat, making it a target for treating metabolic disorders. Here, we investigate how purine nucleotides inhibit respiration uncoupling by UCP1. Our molecular simulations predict that GDP and GTP bind UCP1 in the common substrate binding site in an upright orientation, where the base moiety interacts with conserved residues R92 and E191. We identify a triplet of uncharged residues, F88/I187/W281, forming hydrophobic contacts with nucleotides. In yeast spheroplast respiration assays, both I187A and W281A mutants increase the fatty acid-induced uncoupling activity of UCP1 and partially suppress the inhibition of UCP1 activity by nucleotides. The F88A/I187A/W281A triple mutant is overactivated by fatty acids even at high concentrations of purine nucleotides. In simulations, E191 and W281 interact with purine but not pyrimidine bases. These results provide a molecular understanding of the selective inhibition of UCP1 by purine nucleotides.

Unwrapping NPT Simulations to Calculate Diffusion Coefficients.

In molecular dynamics simulations in the NPT ensemble at constant pressure, the size and shape of the periodic simulation box fluctuate with time. For particle images far from the origin, the rescaling of the box by the barostat results in unbounded position displacements. Special care is thus required when a particle trajectory is unwrapped from a projection into the central box under periodic boundary conditions to a trajectory in full three-dimensional space, e.g., for the calculation of translational diffusion coefficients. Here, we review and compare different schemes in use for trajectory unwrapping. We also specify the corresponding rewrapping schemes to put an unwrapped trajectory back into the central box. On this basis, we then identify a scheme for the calculation of diffusion coefficients from NPT simulations, which is a primary application of trajectory unwrapping. In this scheme, the wrapped and unwrapped trajectory are mutually consistent and their statistical properties are preserved. We conclude with advice on best practice for the consistent unwrapping of constant-pressure simulation trajectories and the calculation of accurate translational diffusion coefficients.

Open-channel structure of a pentameric ligand-gated ion channel reveals a mechanism of leaflet-specific phospholipid modulation.

Pentameric ligand-gated ion channels (pLGICs) mediate synaptic transmission and are sensitive to their lipid environment. The mechanism of phospholipid modulation of any pLGIC is not well understood. We demonstrate that the model pLGIC, ELIC (Erwinia ligand-gated ion channel), is positively modulated by the anionic phospholipid, phosphatidylglycerol, from the outer leaflet of the membrane. To explore the mechanism of phosphatidylglycerol modulation, we determine a structure of ELIC in an open-channel conformation. The structure shows a bound phospholipid in an outer leaflet site, and structural changes in the phospholipid binding site unique to the open-channel. In combination with streamlined alchemical free energy perturbation calculations and functional measurements in asymmetric liposomes, the data support a mechanism by which an anionic phospholipid stabilizes the activated, open-channel state of a pLGIC by specific, state-dependent binding to this site.

Consistent Picture of Phosphate-Divalent Cation Binding from Models with Implicit and Explicit Electronic Polarization.

The binding of divalent cations to the ubiquitous phosphate group is essential for a number of key biological processes, such as DNA compaction, RNA folding, or interactions of some proteins with membranes. Yet, probing their binding sites, modes, and associated binding free energy is a challenge for both experiments and simulations. In simulations, standard force fields strongly overestimate the interaction between phosphate groups and divalent cations. Here, we examine how different strategies to include electronic polarization effects in force fields─implicitly, through the use of scaled charges or pair-specific Lennard-Jones parameters, or explicitly, with the polarizable force fields Drude and AMOEBA─capture the interactions of a model phosphate compound, dimethyl phosphate, with calcium and magnesium divalent cations. We show that both implicit and explicit approaches, when carefully parameterized, are successful in capturing the overall binding free energy and that common trends emerge from the comparison of different simulation approaches. Overall, the binding is very moderate, slightly weaker for Ca2+ than Mg2+, and the solvent-shared ion pair is slightly more stable than the contact monodentate ion pair. The bidentate ion pair is higher in energy (or even fully unstable for Mg2+). Our results thus suggest practical ways to capture the divalent cations with biomolecular phosphate groups in complex biochemical systems. In particular, the computational efficiency of implicit models makes them ideally suited for large-scale simulations of biological assemblies, with improved accuracy compared to state-of-the-art fixed-charge force fields.

Symmetry-Adapted Restraints for Binding Free Energy Calculations.

Binding free energy calculations rely critically on a precise definition of the bound state and well-designed ligand restraints to ensure that binding free energy calculations converge rapidly and yield estimates of well-defined thermodynamic quantities. The distance-to-bound-configuration (DBC) is a single variable that can precisely delineate the bound state of a ligand including translational, rotational and conformational degrees of freedom and has been successfully used to capture binding modes with complex geometries. DBC is defined as the root-mean-square deviation (RMSD) of ligand coordinates in the frame of reference of the binding site. In the special case where the ligand features symmetry-equivalent atoms, a standard RMSD arbitrarily distinguishes equivalent poses, mixing equivalent and nonequivalent degrees of freedom, and preventing the precise delineation of the bound state ensemble, which negates the benefits of defining a flat-bottom binding restraint. To remedy this, we introduce a symmetry-adapted DBC coordinate where the RMSD is minimized over permutations of equivalent ligand atoms. This coordinate is implemented in a portable software library, the Collective Variables Module. We tested the approach by computing the absolute binding free energy of benzene to the engineered site of a mutant lysozyme (L99A/M102H) using alchemical free energy perturbation. We found that the symmetry-adapted restraint leads to well-behaved convergence of both the decoupling free energy in the binding site and the restrained free energy in the gas phase, recovering the affinity computed using a classic center-of-mass restraint. Thus, symmetry-adapted DBC seamlessly generalizes the benefits of DBC restraints to the case of symmetric ligands. The underlying symmetric RMSD coordinate can also be used for analyzing or biasing simulations in other contexts than affinity predictions.

Human Learning for Molecular Simulations: The Collective Variables Dashboard in VMD.

The Collective Variables Dashboard is a software tool for real-time, seamless exploration of molecular structures and trajectories in a customizable space of collective variables. The Dashboard arises from the integration of the Collective Variables Module (also known as Colvars) with the visualization software VMD, augmented with a fully discoverable graphical interface offering interactive workflows for the design and analysis of collective variables. Typical use cases include a priori design of collective variables for enhanced sampling and free energy simulations as well as analysis of any type of simulation or collection of structures in a collective variable space. A combination of those cases commonly occurs when preliminary simulations, biased or unbiased, reveal that an optimized set of collective variables is necessary to improve sampling in further simulations. Then the Dashboard provides an efficient way to intuitively explore the space of likely collective variables, validate them on existing data, and use the resulting collective variable definitions directly in further biased simulations using the Collective Variables Module. Visualization of biasing energies and forces is proposed to help analyze or plan biased simulations. We illustrate the use of the Dashboard on two applications: discovering coordinates to describe ligand unbinding from a protein binding site and designing volume-based variables to bias the hydration of a transmembrane pore.

Multiscale Computational Study of the Conformation of the Full-Length Intrinsically Disordered Protein MeCP2.

The malfunction of the methyl-CpG binding protein 2 (MeCP2) is associated with the Rett syndrome, one of the most common causes of cognitive impairment in females. MeCP2 is an intrinsically disordered protein (IDP), making its experimental characterization a challenge. There is currently no structure available for the full-length MeCP2 in any of the databases, and only the structure of its MBD domain has been solved. We used this structure to build a full-length model of MeCP2 by completing the rest of the protein via ab initio modeling. Using a combination of all-atom and coarse-grained simulations, we characterized its structure and dynamics as well as the conformational space sampled by the ID and transcriptional repression domain (TRD) domains in the absence of the rest of the protein. The present work is the first computational study of the full-length protein. Two main conformations were sampled in the coarse-grained simulations: a globular structure similar to the one observed in the all-atom force field and a two-globule conformation. Our all-atom model is in good agreement with the available experimental data, predicting amino acid W104 to be buried, amino acids R111 and R133 to be solvent-accessible, and having a 4.1% α-helix content, compared to the 4% found experimentally. Finally, we compared the model predicted by AlphaFold to our Modeller model. The model was not stable in water and underwent further folding. Together, these simulations provide a detailed (if perhaps incomplete) conformational ensemble of the full-length MeCP2, which is compatible with experimental data and can be the basis of further studies, e.g., on mutants of the protein or its interactions with its biological partners.

Fast and Accurate Multidimensional Free Energy Integration.

Enhanced sampling and free energy calculation algorithms of the thermodynamic integration family (such as the adaptive biasing force (ABF) method) are not based on the direct computation of a free energy surface but rather of its gradient. Integrating the free energy surface is nontrivial in dimensions higher than one. Here, the author introduces a flexible, portable implementation of a Poisson equation formalism to integrate free energy surfaces from estimated gradients in dimensions 2 and 3 using any combination of periodic and nonperiodic (Neumann) boundary conditions. The algorithm is implemented in portable C++ and provided as a standalone tool that can be used to integrate multidimensional gradient fields estimated on a grid using any algorithm, such as umbrella integration as a post-treatment of umbrella sampling simulations. It is also included in the implementation of ABF (and its extended-system variant eABF) in the Collective Variables Module, enabling the seamless computation of multidimensional free energy surfaces within ABF and eABF simulations. A Python-based analysis toolchain is provided to easily plot and analyze multidimensional ABF simulation results, including metrics to assess their convergence. The Poisson integration algorithm can also be used to perform Helmholtz decomposition of noisy gradient estimates on the fly, resulting in an efficient implementation of the projected ABF (pABF) method proposed by Leliévre and co-workers. In numerical tests, pABF is found to lead to faster convergence with respect to ABF in simple cases of low intrinsic dimension but seems detrimental to convergence in a more realistic case involving degenerate coordinates and hidden barriers due to slower exploration. This suggests that variance reduction schemes do not always yield convergence improvements when applied to enhanced sampling methods.

Mechanistic Insights on Heme-to-Heme Transmembrane Electron Transfer Within NADPH Oxydases From Atomistic Simulations.

NOX5 is a member of the NADPH oxidase family which is dedicated to the production of reactive oxygen species. The molecular mechanisms governing transmembrane electron transfer (ET) that permits to shuttle electrons over the biological membrane have remained elusive for a long time. Using computer simulations, we report conformational dynamics of NOX5 embedded within a realistic membrane environment. We assess the stability of the protein within the membrane and monitor the existence of cavities that could accommodate dioxygen molecules. We investigate the heme-to-heme electron transfer. We find a reaction free energy of a few tenths of eV (ca. -0.3 eV) and a reorganization free energy of around 1.1 eV (0.8 eV after including electrostatic induction corrections). The former indicates thermodynamically favorable ET, while the latter falls in the expected values for transmembrane inter-heme ET. We estimate the electronic coupling to fall in the range of the µeV. We identify electron tunneling pathways showing that not only the W378 residue is playing a central role, but also F348. Finally, we reveal the existence of two connected O2-binding pockets near the outer heme with fast exchange between the two sites on the nanosecond timescale. We show that when the terminal heme is reduced, O2 binds closer to it, affording a more efficient tunneling pathway than when the terminal heme is oxidized, thereby providing an efficient mechanism to catalyze superoxide production in the final step. Overall, our study reveals some key molecular mechanisms permitting reactive oxygen species production by NOX5 and paves the road for further investigation of ET processes in the wide family of NADPH oxidases by computer simulations.

Binding of divalent cations to acetate: molecular simulations guided by Raman spectroscopy.

In spite of the biological importance of the binding of Zn2+, Ca2+, and Mg2+ to the carboxylate group, cation-acetate binding affinities and binding modes remain actively debated. Here, we report the first use of Raman multivariate curve resolution (Raman-MCR) vibrational spectroscopy to obtain self-consistent free and bound metal acetate spectra and one-to-one binding constants, without the need to invoke any a priori assumptions regarding the shapes of the corresponding vibrational bands. The experimental results, combined with classical molecular dynamics simulations with a force field effectively accounting for electronic polarization via charge scaling and ab initio simulations, indicate that the measured binding constants pertain to direct (as opposed to water separated) ion pairing. The resulting binding constants do not scale with cation size, as the binding constant to Zn2+ is significantly larger than that to either Mg2+ or Ca2+, although Zn2+ and Mg2+ have similar radii that are about 25% smaller than Ca2+. Remaining uncertainties in the metal acetate binding free energies are linked to fundamental ambiguities associated with identifying the range of structures pertaining to non-covalently bound species.

Scalable molecular dynamics on CPU and GPU architectures with NAMD.

NAMDis a molecular dynamics program designed for high-performance simulations of very large biological objects on CPU- and GPU-based architectures. NAMD offers scalable performance on petascale parallel supercomputers consisting of hundreds of thousands of cores, as well as on inexpensive commodity clusters commonly found in academic environments. It is written in C++ and leans on Charm++ parallel objects for optimal performance on low-latency architectures. NAMD is a versatile, multipurpose code that gathers state-of-the-art algorithms to carry out simulations in apt thermodynamic ensembles, using the widely popular CHARMM, AMBER, OPLS, and GROMOS biomolecular force fields. Here, we review the main features of NAMD that allow both equilibrium and enhanced-sampling molecular dynamics simulations with numerical efficiency. We describe the underlying concepts utilized by NAMD and their implementation, most notably for handling long-range electrostatics; controlling the temperature, pressure, and pH; applying external potentials on tailored grids; leveraging massively parallel resources in multiple-copy simulations; and hybrid quantum-mechanical/molecular-mechanical descriptions. We detail the variety of options offered by NAMD for enhanced-sampling simulations aimed at determining free-energy differences of either alchemical or geometrical transformations and outline their applicability to specific problems. Last, we discuss the roadmap for the development of NAMD and our current efforts toward achieving optimal performance on GPU-based architectures, for pushing back the limitations that have prevented biologically realistic billion-atom objects to be fruitfully simulated, and for making large-scale simulations less expensive and easier to set up, run, and analyze. NAMD is distributed free of charge with its source code at www.ks.uiuc.edu.

A Streamlined, General Approach for Computing Ligand Binding Free Energies and Its Application to GPCR-Bound Cholesterol.

The theory of receptor-ligand binding equilibria has long been well-established in biochemistry, and was primarily constructed to describe dilute aqueous solutions. Accordingly, few computational approaches have been developed for making quantitative predictions of binding probabilities in environments other than dilute isotropic solution. Existing techniques, ranging from simple automated docking procedures to sophisticated thermodynamics-based methods, have been developed with soluble proteins in mind. Biologically and pharmacologically relevant protein-ligand interactions often occur in complex environments, including lamellar phases like membranes and crowded, nondilute solutions. Here, we revisit the theoretical bases of ligand binding equilibria, avoiding overly specific assumptions that are nearly always made when describing receptor-ligand binding. Building on this formalism, we extend the asymptotically exact Alchemical Free Energy Perturbation technique to quantifying occupancies of sites on proteins in a complex bulk, including phase-separated, anisotropic, or nondilute solutions, using a thermodynamically consistent and easily generalized approach that resolves several ambiguities of current frameworks. To incorporate the complex bulk without overcomplicating the overall thermodynamic cycle, we simplify the common approach for ligand restraints by using a single distance-from-bound-configuration (DBC) ligand restraint during AFEP decoupling from protein. DBC restraints should be generalizable to binding modes of most small molecules, even those with strong orientational dependence. We apply this approach to compute the likelihood that membrane cholesterol binds to known crystallographic sites on three GPCRs (ß2-adrenergic, 5HT-2B, and µ-opioid) at a range of concentrations. Nonideality of cholesterol in a binary cholesterol:phosphatidylcholine (POPC) bilayer is characterized and consistently incorporated into the interpretation. We find that the three sites exhibit very different affinities for cholesterol: The site on the adrenergic receptor is predicted to be high affinity, with 50% occupancy for 1:109 CHOL:POPC mixtures. The sites on the 5HT-2B and µ-opioid receptor are predicted to be lower affinity, with 50% occupancy for 1:103 CHOL:POPC and 1:102 CHOL:POPC, respectively. These results could not have been predicted from the crystal structures alone.

Molecular Mechanics Parameterization of Anesthetic Molecules.

Anesthetic drug molecules are being increasingly studied through the use of computational methods such as molecular dynamics (MD). Molecular mechanics force fields require the investigator to supply parameters for the force field equation, which are not available for novel molecules. Careful selection of these parameters is critical for simulations to produce meaningful results. Therefore, this chapter presents a state-of-the-art method for determining these parameters by comparison to quantum mechanics calculations and experimental quantities. Ketamine is used as an example to demonstrate the process.

New Coarse Variables for the Accurate Determination of Standard Binding Free Energies.

To improve sampling of the configurational entropy change upon protein-ligand binding, we have introduced a new set of coarse variables describing the relative orientation and position of the ligand via a global macromolecular orientational procedure, onto which geometrical restraints are applied. Evaluating the potential of mean force for the different coarse variables, the experimental standard binding free energy for three decapeptides associated with the SH3 domain of the Abl kinase is reproduced quantitatively.

A membrane-inserted structural model of the yeast mitofusin Fzo1.

Mitofusins are large transmembrane GTPases of the dynamin-related protein family, and are required for the tethering and fusion of mitochondrial outer membranes. Their full-length structures remain unknown, which is a limiting factor in the study of outer membrane fusion. We investigated the structure and dynamics of the yeast mitofusin Fzo1 through a hybrid computational and experimental approach, combining molecular modelling and all-atom molecular dynamics simulations in a lipid bilayer with site-directed mutagenesis and in vivo functional assays. The predicted architecture of Fzo1 improves upon the current domain annotation, with a precise description of the helical spans linked by flexible hinges, which are likely of functional significance. In vivo site-directed mutagenesis validates salient aspects of this model, notably, the long-distance contacts and residues participating in hinges. GDP is predicted to interact with Fzo1 through the G1 and G4 motifs of the GTPase domain. The model reveals structural determinants critical for protein function, including regions that may be involved in GTPase domain-dependent rearrangements.

Smoothed Biasing Forces Yield Unbiased Free Energies with the Extended-System Adaptive Biasing Force Method.

We report a theoretical description and numerical tests of the extended-system adaptive biasing force method (eABF), together with an unbiased estimator of the free energy surface from eABF dynamics. Whereas the original ABF approach uses its running estimate of the free energy gradient as the adaptive biasing force, eABF is built on the idea that the exact free energy gradient is not necessary for efficient exploration, and that it is still possible to recover the exact free energy separately with an appropriate estimator. eABF does not directly bias the collective coordinates of interest, but rather fictitious variables that are harmonically coupled to them; therefore is does not require second derivative estimates, making it easily applicable to a wider range of problems than ABF. Furthermore, the extended variables present a smoother, coarse-grain-like sampling problem on a mollified free energy surface, leading to faster exploration and convergence. We also introduce CZAR, a simple, unbiased free energy estimator from eABF trajectories. eABF/CZAR converges to the physical free energy surface faster than standard ABF for a wide range of parameters.

A Novel Bifunctional Alkylphenol Anesthetic Allows Characterization of γ-Aminobutyric Acid, Type A (GABAA), Receptor Subunit Binding Selectivity in Synaptosomes.

Propofol, an intravenous anesthetic, is a positive modulator of the GABAA receptor, but the mechanistic details, including the relevant binding sites and alternative targets, remain disputed. Here we undertook an in-depth study of alkylphenol-based anesthetic binding to synaptic membranes. We designed, synthesized, and characterized a chemically active alkylphenol anesthetic (2-((prop-2-yn-1-yloxy)methyl)-5-(3-(trifluoromethyl)-3H-diazirin-3-yl)phenol, AziPm-click (1)), for affinity-based protein profiling (ABPP) of propofol-binding proteins in their native state within mouse synaptosomes. The ABPP strategy captured ∼4% of the synaptosomal proteome, including the unbiased capture of five α or ß GABAA receptor subunits. Lack of γ2 subunit capture was not due to low abundance. Consistent with this, independent molecular dynamics simulations with alchemical free energy perturbation calculations predicted selective propofol binding to interfacial sites, with higher affinities for α/ß than γ-containing interfaces. The simulations indicated hydrogen bonding is a key component leading to propofol-selective binding within GABAA receptor subunit interfaces, with stable hydrogen bonds observed between propofol and α/ß cavity residues but not γ cavity residues. We confirmed this by introducing a hydrogen bond-null propofol analogue as a protecting ligand for targeted-ABPP and observed a lack of GABAA receptor subunit protection. This investigation demonstrates striking interfacial GABAA receptor subunit selectivity in the native milieu, suggesting that asymmetric occupancy of heteropentameric ion channels by alkylphenol-based anesthetics is sufficient to induce modulation of activity.

Evolution of Pentameric Ligand-Gated Ion Channels: Pro-Loop Receptors.

Pentameric ligand-gated ion channels (pLGICs) are ubiquitous neurotransmitter receptors in Bilateria, with a small number of known prokaryotic homologues. Here we describe a new inventory and phylogenetic analysis of pLGIC genes across all kingdoms of life. Our main finding is a set of pLGIC genes in unicellular eukaryotes, some of which are metazoan-like Cys-loop receptors, and others devoid of Cys-loop cysteines, like their prokaryotic relatives. A number of such "Cys-less" receptors also appears in invertebrate metazoans. Together, those findings draw a new distribution of pLGICs in eukaryotes. A broader distribution of prokaryotic channels also emerges, including a major new archaeal taxon, Thaumarchaeota. More generally, pLGICs now appear nearly ubiquitous in major taxonomic groups except multicellular plants and fungi. However, pLGICs are sparsely present in unicellular taxa, suggesting a high rate of gene loss and a non-essential character, contrasting with their essential role as synaptic receptors of the bilaterian nervous system. Multiple alignments of these highly divergent sequences reveal a small number of conserved residues clustered at the interface between the extracellular and transmembrane domains. Only the "Cys-loop" proline is absolutely conserved, suggesting the more fitting name "Pro loop" for that motif, and "Pro-loop receptors" for the superfamily. The infered molecular phylogeny shows a Cys-loop and a Cys-less clade in eukaryotes, both containing metazoans and unicellular members. This suggests new hypotheses on the evolutionary history of the superfamily, such as a possible origin of the Cys-loop cysteines in an ancient unicellular eukaryote. Deeper phylogenetic relationships remain uncertain, particularly around the split between bacteria, archaea, and eukaryotes.