A millisecond coarse-grained simulation approach to decipher allosteric cannabinoid binding at the glycine receptor α1

Glycine receptors (GlyR) are regulated by small-molecule binding at several allosteric sites. Cannabinoids like tetrahydrocannabinol (THC) and N-arachidonyl-ethanol-amide (AEA) potentiate the GlyR response but their mechanism of action is not fully established. By combining millisecond coarse-grained (CG) MD simulations powered by Martini 3 with backmapping to all-atom representations, we have characterized the cannabinoid-binding site(s) at the zebrafish GlyR-α1 active state with atomic resolution. Based on hundreds of thousand ligand-binding events, we find that cannabinoids bind to the transmembrane domain of the receptor at both intrasubunit and intersubunit sites. For THC, the intrasubunit binding mode predicted in simulation is in excellent agreement with recent cryo-EM structures, while intersubunit binding recapitulates in full previous mutagenesis experiments. Intriguingly, AEA is predicted to bind at the same intersubunit site despite the strikingly different chemistry. Statistical analyses of the ligand-receptor interactions highlight potentially relevant residues for GlyR potentiation, offering experimentally testable predictions. The predictions for AEA have been validated by electrophysiology recordings of rationally designed mutants. The results highlight the existence of multiple cannabinoid-binding sites for the allosteric regulation of GlyR and put forward an effective strategy for the identification and structural characterization of allosteric binding sites.

Modeling Chromatography Binding through Molecular Dynamics Simulations with Resin Fragments

Accurate atomistic modeling of the interactions of a chromatography resin with a solute can inform the selection of purification conditions for a product, an important problem in the biotech and pharmaceutical industries. We present a molecular dynamics simulation-based approach for the qualitative prediction of interaction sites (specificity) and retention times (affinity) of a protein for a given chromatography resin. We mimicked the resin with an unrestrained ligand composed of the resin headgroup coupled with successively larger fragments of the agarose backbone. The interactions of the ligand with the protein are simulated in an explicit solvent using the Replica Exchange Molecular Dynamics enhanced sampling approach in conjunction with Hydrogen Mass Repartitioning (REMD-HMR). We computed the ligand interaction surface from the simulation trajectories and correlated the features of the interaction surface with experimentally determined retention times. The simulation and analysis protocol were first applied to a series of ubiquitin mutants for which retention times on Capto MMC resin are available. The ubiquitin simulations helped identify the optimal ligand that was used in subsequent simulations on six proteins for which Capto MMC elution times are available. For each of the six proteins, we computed the interaction surface and characterized it in terms of a range of simulation-averaged residue-level physicochemical descriptors. Modeling of the salt concentrations required for elution with respect to the descriptors resulted in a linear fit in terms of aromaphilicity and Kyte-Doolittle hydrophobicity that was robust to outliers, showed high correlation, and correctly ranked the protein elution order. The physics-based model building approach described here does not require a large experimental data set and can be readily applied to different resins and diverse biomolecules.

Breakdown of Langmuir Adsorption Isotherm in Small Closed Systems

For more than a century, monolayer adsorptions in which adsorbate molecules and adsorbing sites behave ideally have been successfully described by Langmuir's adsorption isotherm. For example, the amount of adsorbed material, as a function of concentration of the material which is not adsorbed, obeys Langmuir's equation. In this paper, we argue that this relation is valid only for macroscopic systems. However, when particle numbers of adsorbate molecules and/or adsorbing sites are small, Langmuir's model fails to describe the chemical equilibrium of the system. This is because the kinetics of forming, or the probability of observing, occupied sites arises from two-body interactions, and as such, ought to include cross-correlations between particle numbers of the adsorbate and adsorbing sites. The effect of these correlations, as reflected by deviations in predicting composition when correlations are ignored, increases with decreasing particle numbers and becomes substantial when only few adsorbate molecules, or adsorbing sites, are present in the system. In addition, any change that augments the fraction of occupied sites at equilibrium (e.g., smaller volume, lower temperature, or stronger adsorption energy) further increases the discrepancy between observed properties of small systems and those predicted by Langmuir's theory. In contrast, for large systems, these cross-correlations become negligible, and therefore when expressing properties involving two-body processes, it is possible to consider independently the concentration of each component. By applying statistical mechanics concepts, we derive a general expression of the equilibrium constant for adsorption. It is also demonstrated that in ensembles in which total numbers of particles are fixed, the magnitudes of fluctuations in particle numbers alone can predict the average chemical composition of the system. Moreover, an alternative adsorption equation, predicting the average fraction of occupied sites from the value of the equilibrium constant, is proposed. All derived relations were tested against results obtained by Monte Carlo simulations.

Hexagonal Lattices of HIV Capsid Proteins Explored by Simulations Based on a Thermodynamically Consistent Model

HIV capsid proteins (CAs) may self-assemble into a variety of shapes under in vivo and in vitro conditions. Here, we employed simulations based on a residue-level coarse-grained (CG) model with full conformational flexibility to investigate hexagonal lattices, which are the underlying structural pattern for CA aggregations. Facilitated by enhanced sampling simulations to rigorously calculate CA dimerization and polymerization affinities, we calibrated our model to reproduce the experimentally measured affinities. Using the calibrated model, we performed unbiased simulations on several large systems consisting of 1512 CA subunits, allowing reversible binding and unbinding of the CAs in a thermodynamically consistent manner. In one simulation, a preassembled hexagonal CA sheet developed spontaneous curvatures reminiscent of those observed in experiments, and the edges of the sheet exhibited local curvatures larger than those of the interior. In other simulations starting with randomly distributed CAs at different concentrations, existing CA assemblies grew by binding free capsomeres to the edges and by merging with other assemblies. At high CA concentrations, rapid establishment of predominant aggregates was followed by much slower adjustments toward more regular hexagonal lattices, with increasing numbers of intact CA hexamers and pentamers being formed. Our approach of adapting a general CG model to specific systems by using experimental binding data represents a practical and effective strategy for simulating and elucidating intricate protein aggregations.

A computational toolbox for the assembly yield of complex and heterogeneous structures

The self-assembly of complex structures from a set of non-identical building blocks is a hallmark of soft matter and biological systems, including protein complexes, colloidal clusters, and DNA-based assemblies. Predicting the dependence of the equilibrium assembly yield on the concentrations and interaction energies of building blocks is highly challenging, owing to the difficulty of computing the entropic contributions to the free energy of the many structures that compete with the ground state configuration. While these calculations yield well known results for spherically symmetric building blocks, they do not hold when the building blocks have internal rotational degrees of freedom. Here we present an approach for solving this problem that works with arbitrary building blocks, including proteins with known structure and complex colloidal building blocks. Our algorithm combines classical statistical mechanics with recently developed computational tools for automatic differentiation. Automatic differentiation allows efficient evaluation of equilibrium averages over configurations that would otherwise be intractable. We demonstrate the validity of our framework by comparison to molecular dynamics simulations of simple examples, and apply it to calculate the yield curves for known protein complexes and for the assembly of colloidal shells.

Dynamics of organophosphate guest encapsulation in heteroleptic coordination cages

Heteroleptic coordination cages allow the design of different host structures that can bind guest molecules within their cavities. In a previous work, the energetics of organophosphate encapsulation in palladium(II)-based heteroleptic coordination cages that differ in terms of their ability to form hydrogen bonds have been investigated [Platzek et al., Endohedrally Functionalized Heteroleptic Coordination Cages for Phosphate Ester Binding, Angew. Chem., Int. Ed. 2022, 61, e2022093]. The present work focuses on the dynamics of this system. Dynamic information is obtained through the application of a Markov state model (MSM) to unbiased multi-microsecond atomistic molecular dynamics simulations of guest binding and release. The MSM reveals that both the bound state and the binding/unbinding pathways are highly dynamic, with different types of interactions mediating the binding of the diphenylphosphate guest. Thus, the simulations highlight the dynamic nature of the nanoconfinement in the host-guest systems, with possible implications for the use of such coordination cages as catalysts.

Insights into the coordination chemistry of antineoplastic doxorubicin with 3d-transition metal ions Zn2+, Cu2+, and VO2+: a study using well-calibrated thermodynamic cycles and chemical interaction quantum chemistry models

We present a computational strategy based on thermodynamic cycles to predict and describe the chemical equilibrium between the 3d-transition metal ions Zn²⁺, Cu²⁺, and VO²⁺ and the widely used antineoplastic drug doxorubicin. Our method involves benchmarking a theoretical protocol to compute gas-phase quantities using DLPNO Coupled-Cluster calculations as reference, followed by estimating solvation contributions to the reaction Gibbs free energies using both explicit partial (micro)solvation steps for charged solutes and neutral coordination complexes, as well as a continuum solvation procedure for all solutes involved in the complexation process. We rationalized the stability of these doxorubicin-metal complexes by inspecting quantities obtained from the topology of their electron densities, particularly the bond critical points and non-covalent interaction index. Our approach allowed us to identify representative species in solution phase, infer the most likely complexation process for each case, and identify key intramolecular interactions involved in the stability of these compounds. To the best of our knowledge, this is the first study reporting thermodynamic constants for the complexation of doxorubicin with transition metal ions. Unlike other methods, our procedure is computationally affordable for medium-sized systems and provides valuable insights even with limited experimental data. Furthermore, it can be extended to describe the complexation process between 3d-transition metal ions and other bioactive ligands.

Endohedrally Functionalized Heteroleptic Coordination Cages for Phosphate Ester Binding

Metallosupramolecular hosts of nanoscopic dimensions, which are able to serve as selective receptors and catalysts, are usually composed of only one type of organic ligand, restricting diversity in terms of cavity shape and functional group decoration. We report a series of heteroleptic [Pd2 A2 B2 ] coordination cages that self-assemble from a library of shape complementary bis-monodentate ligands in a non-statistical fashion. Ligands A feature an inward pointing NH function, able to engage in hydrogen bonding and amenable to being functionalized with amide and alkyl substituents. Ligands B comprise tricyclic aromatic backbones of different shape and electronic situation. The obtained heteroleptic coordination cages were investigated for their ability to bind phosphate diesters as guests. All-atom molecular dynamics (MD) simulations in explicit solvent were conducted to understand the mechanistic relationships behind the experimentally determined guest affinities.

Small molecules targeting the disordered transactivation domain of the androgen receptor induce the formation of collapsed helical states

Intrinsically disordered proteins, which do not adopt well-defined structures under physiological conditions, are implicated in many human diseases. Small molecules that target the disordered transactivation domain of the androgen receptor have entered human trials for the treatment of castration-resistant prostate cancer (CRPC), but no structural or mechanistic rationale exists to explain their inhibition mechanisms or relative potencies. Here, we utilize all-atom molecular dynamics computer simulations to elucidate atomically detailed binding mechanisms of the compounds EPI-002 and EPI-7170 to the androgen receptor. Our simulations reveal that both compounds bind at the interface of two transiently helical regions and induce the formation of partially folded collapsed helical states. We find that EPI-7170 binds androgen receptor more tightly than EPI-002 and we identify a network of intermolecular interactions that drives higher affinity binding. Our results suggest strategies for developing more potent androgen receptor inhibitors and general strategies for disordered protein drug design.

Discrete stochastic models of SELEX: Aptamer capture probabilities and protocol optimization

Antibodies are important biomolecules that are often designed to recognize target antigens. However, they are expensive to produce and their relatively large size prevents their transport across lipid membranes. An alternative to antibodies is aptamers, short ([Formula: see text] bp) oligonucleotides (and amino acid sequences) with specific secondary and tertiary structures that govern their affinity to specific target molecules. Aptamers are typically generated via solid phase oligonucleotide synthesis before selection and amplification through Systematic Evolution of Ligands by EXponential enrichment (SELEX), a process based on competitive binding that enriches the population of certain strands while removing unwanted sequences, yielding aptamers with high specificity and affinity to a target molecule. Mathematical analyses of SELEX have been formulated in the mass action limit, which assumes large system sizes and/or high aptamer and target molecule concentrations. In this paper, we develop a fully discrete stochastic model of SELEX. While converging to a mass-action model in the large system-size limit, our stochastic model allows us to study statistical quantities when the system size is small, such as the probability of losing the best-binding aptamer during each round of selection. Specifically, we find that optimal SELEX protocols in the stochastic model differ from those predicted by a deterministic model.

Ion-bridges and lipids drive aggregation of same-charge nanoparticles on lipid membranes

The control of the aggregation of biomedical nanoparticles (NP) in physiological conditions is crucial as clustering may change completely the way they interact with the biological environment. Here we show that Au nanoparticles, functionalized by an anionic, amphiphilic shell, spontaneously aggregate in fluid zwitterionic lipid bilayers. We use molecular dynamics and enhanced sampling techniques to disentangle the short-range and long-range driving forces of aggregation. At short inter-particle distances, ion-mediated, charge-charge interactions (ion bridging) stabilize the formation of large NP aggregates, as confirmed by cryo-electron microscopy. Lipid depletion and membrane curvature are the main membrane deformations driving long-range NP-NP attraction. Ion bridging, lipid depletion, and membrane curvature stem from the configurational flexibility of the nanoparticle shell. Our simulations show, more in general, that the aggregation of same-charge membrane inclusions can be expected as a result of intrinsically nanoscale effects taking place at the NP-NP and NP-bilayer soft interfaces.

Binding reactions at finite systems

A perpetual yearn exists among computational scientists to scale down the size of physical systems, a desire shared as well with experimentalists able to track single molecules. A question then arises whether averages observed at small systems are the same as those observed at large or macroscopic systems. Utilizing statistical-mechanics formulations in ensembles in which the total numbers of particles are fixed, we demonstrate that properties of binding reactions are not homogeneous functions. This means that averages of intensive parameters, such as the concentration of the bound-state, at finite systems are different than those at large systems. The discrepancy increases with decreasing temperature, volume, and to some extent, numbers of particles. As perplexing as it may sound, despite variations in average quantities, extracting the equilibrium constant from systems of different sizes does yield the same value. The reason is that correlations in reactants' concentrations ought to be accounted for in the expression of the equilibrium constant, being negligible at large-scale but significant at small-scale. Similar arguments pertain to the calculations of the reaction rate constants, more specifically, the bimolecular rate of the forward reaction is related to the average of the product (and not to the product of the averages) of the reactants' concentrations. Furthermore, we derive relations aiming to predict the composition only from the equilibrium constant and the system's size. All predictions are validated by Monte-Carlo and molecular dynamics simulations. An important consequence of these findings is that the expression of the equilibrium constant at finite systems is not dictated solely by the chemical equation of the reaction but requires knowledge of the elementary processes involved.

Computational Tools for Accurate Binding Free-Energy Prediction

A quantitative thermodynamic understanding of the noncovalent association of (bio)molecules is of central importance in molecular life sciences. An important quantity characterizing (bio)molecular association is the binding affinity or absolute binding free energy. In recent years, the computational prediction of absolute binding free energies has evolved considerably in terms of accuracy, computational speed, and user-friendliness. In this chapter, we first give an overview of how absolute free energies are defined and how they can be determined with computational means. We proceed with an outline of the theoretical basis of the two most reliable methods, potential of mean force, and double decoupling calculations. In particular, we describe how the sampling problem can be alleviated by application of restraints. Finally, we provide step-by-step instructions of how to set up corresponding molecular simulations with a commonly employed molecular dynamics simulation engine.

Evolution of Thyroglobulin Loop Kinetics in EpCAM

Epithelial cell-activating molecule (EpCAM) is an important cancer biomarker and therapeutic target given its elevated expression in epithelial cancers. EpCAM is a type I transmembrane protein that forms cis-dimers along the thyroglobulin type-1A-like domain (TYD) in the extracellular region. The thyroglobulin loop (TY loop) within the TYD is structurally dynamic in the monomer state of human EpCAM, binding reversibly to a TYD site. However, it is not known if this flexibility is prevalent across different species. Here, we conduct over 17 μs of all-atom molecular dynamics simulations to study EpCAM TY loop kinetics of five different species, including human, mouse, chicken, frog, and fish. We find that the TY loop remains dynamic across evolution. In addition to the TYD binding site, we discover a second binding site for the TY loop in the C-terminal domain (CTD). Calculations of the dissociation rate constants from the simulation trajectories suggest a differential binding pattern of fish EpCAM and other organisms. Whereas fish TY loop has comparable binding for both TYD and CTD sites, the TY loops of other species preferably bind the TYD site. A hybrid construct of fish EpCAM with human TY loop restores the TYD binding preference, suggesting robust effects of the TY loop sequence on its dynamic behavior. Our findings provide insights into the structural dynamics of EpCAM and its implication in physiological functions.

Design of peptides with strong binding affinity to poly(methyl methacrylate) resin by use of molecular simulation-based materials informatics

Peptides with strong binding affinities for poly(methyl methacrylate) (PMMA) resin were designed by use of materials informatics technology based on molecular dynamics simulation for the purpose of covering the resin surface with adhesive peptides, which were expected to result in eco-friendly and biocompatible biomaterials. From the results of binding affinity obtained with this molecular simulation, it was confirmed that experimental values could be predicted with errors <10%. By analyzing the simulation data with the response-surface method, we found that three peptides (RWWRPWW, EWWRPWR, and RWWRPWR), which consist of arginine (R), tryptophan (W), and proline (P), have strong binding affinity to the PMMA resin. These amino acids were effective because arginine and tryptophan have strong binding affinities for methoxycarbonyl groups and methyl groups, which are the main constituents of the PMMA resin, and proline stabilizes the flat zigzag structures of the peptides in water. The strong binding affinities of the three peptides were confirmed by experiments (surface plasmon resonance methods).

Physics of the Nuclear Pore Complex: Theory, Modeling and Experiment

The hallmark of eukaryotic cells is the nucleus that contains the genome, enclosed by a physical barrier known as the nuclear envelope (NE). On the one hand, this compartmentalization endows the eukaryotic cells with high regulatory complexity and flexibility. On the other hand, it poses a tremendous logistic and energetic problem of transporting millions of molecules per second across the nuclear envelope, to facilitate their biological function in all compartments of the cell. Therefore, eukaryotes have evolved a molecular "nanomachine" known as the Nuclear Pore Complex (NPC). Embedded in the nuclear envelope, NPCs control and regulate all the bi-directional transport between the cell nucleus and the cytoplasm. NPCs combine high molecular specificity of transport with high throughput and speed, and are highly robust with respect to molecular noise and structural perturbations. Remarkably, the functional mechanisms of NPC transport are highly conserved among eukaryotes, from yeast to humans, despite significant differences in the molecular components among various species. The NPC is the largest macromolecular complex in the cell. Yet, despite its significant complexity, it has become clear that its principles of operation can be largely understood based on fundamental physical concepts, as have emerged from a combination of experimental methods of molecular cell biology, biophysics, nanoscience and theoretical and computational modeling. Indeed, many aspects of NPC function can be recapitulated in artificial mimics with a drastically reduced complexity compared to biological pores. We review the current physical understanding of the NPC architecture and function, with the focus on the critical analysis of experimental studies in cells and artificial NPC mimics through the lens of theoretical and computational models. We also discuss the connections between the emerging concepts of NPC operation and other areas of biophysics and bionanotechnology.

Building intuition for binding free energy calculations: Bound state definition, restraints, and symmetry

The theory behind computation of absolute binding free energies using explicit-solvent molecular simulations is well-established, yet somewhat complex, with counter-intuitive aspects. This leads to frequent frustration, common misconceptions, and sometimes erroneous numerical treatment. To improve this, we present the main practically relevant segments of the theory with constant reference to physical intuition. We pinpoint the role of the implicit or explicit definition of the bound state (or the binding site) to make a robust link between an experimental measurement and a computational result. We clarify the role of symmetry and discuss cases where symmetry number corrections have been misinterpreted. In particular, we argue that symmetry corrections as classically presented are a source of confusion and could be advantageously replaced by restraint free energy contributions. We establish that contrary to a common intuition, partial or missing sampling of some modes of symmetric bound states does not affect the calculated decoupling free energies. Finally, we review these questions and pitfalls in the context of a few common practical situations: binding to a symmetric receptor (equivalent binding sites), binding of a symmetric ligand (equivalent poses), and formation of a symmetric complex, in the case of homodimerization.

Modeling protein association from homogeneous to mixed environments: A reaction-diffusion dynamics approach

Protein-protein association in vivo occur in a crowded and complex environment. Theoretical models based on hard-core repulsion predict stabilization of the product under crowded conditions. Soft interactions, on the contrary, can either stabilize or destabilize the product formation. Here we modeled protein association in presence of crowders of varying size, shape, interaction potential and used different mixing parameters for constituent crowders to study the influence on the association reaction. It was found that size is a more dominant factor in crowder-induced stabilization than the shape. Furthermore, in a mixture of crowders having different sizes but identical interaction potential, the change of free energy is additive of the free energy changes produced by individual crowders. However, the free energy change is not additive if two crowders of same size interact via different interaction potentials. These findings provide a systematic understanding of crowding influences in heterogeneous medium.

Physical modeling of multivalent interactions in the nuclear pore complex

In the nuclear pore complex, intrinsically disordered proteins (FG Nups), along with their interactions with more globular proteins called nuclear transport receptors (NTRs), are vital to the selectivity of transport into and out of the cell nucleus. Although such interactions can be modeled at different levels of coarse graining, in vitro experimental data have been quantitatively described by minimal models that describe FG Nups as cohesive homogeneous polymers and NTRs as uniformly cohesive spheres, in which the heterogeneous effects have been smeared out. By definition, these minimal models do not account for the explicit heterogeneities in FG Nup sequences, essentially a string of cohesive and noncohesive polymer units, and at the NTR surface. Here, we develop computational and analytical models that do take into account such heterogeneity in a minimal fashion and compare them with experimental data on single-molecule interactions between FG Nups and NTRs. Overall, we find that the heterogeneous nature of FG Nups and NTRs does play a role in determining equilibrium binding properties but is of much greater significance when it comes to unbinding and binding kinetics. Using our models, we predict how binding equilibria and kinetics depend on the distribution of cohesive blocks in the FG Nup sequences and of the binding pockets at the NTR surface, with multivalency playing a key role. Finally, we observe that single-molecule binding kinetics has a rather minor influence on the diffusion of NTRs in polymer melts consisting of FG-Nup-like sequences.

Statistical physics of DNA hybridization

Deoxyribonucleic acid (DNA) hybridization is at the heart of countless biological and biotechnological processes. Its theoretical modeling played a crucial role, since it has enabled extracting the relevant thermodynamic parameters from systematic measurements of DNA melting curves. In this article, we propose a framework based on statistical physics to describe DNA hybridization and melting in an arbitrary mixture of DNA strands. In particular, we are able to analytically derive closed expressions of the system partition functions for any number N of strings and explicitly calculate them in two paradigmatic situations: (i) a system made of self-complementary sequences and (ii) a system comprising two mutually complementary sequences. We derive the melting curve in the thermodynamic limit (N→∞) of our description, which provides a full justification for the extra entropic contribution that in classic hybridization modeling was required to correctly describe within the same framework the melting of sequences either self-complementary or not. We thus provide a thorough study comprising limit cases and alternative approaches showing how our framework can give a comprehensive view of hybridization and melting phenomena.

Thermodynamic driving forces of guest confinement in a photoswitchable cage

Photoswitchable cages that confine small guest molecules inside their cavities offer a way to control the binding/unbinding process through irradiation with light of different wavelengths. However, detailed characterization of the structural and thermodynamic consequences of photoswitching is very challenging to achieve by experiments alone. Thus, all-atom molecular dynamics (MD) simulations were carried out to gain insight into the relationship between the structure and binding affinity. Binding free energies of the B₁₂F₁₂^2- guest were obtained for all photochemically accessible forms of a photoswitchable dithienylethene (DTE) based coordination cage. The MD simulations show that successive photo-induced closure of the four individual DTE ligands that form the cage gradually decreases the binding affinity. Closure of the first ligand significantly lowers the unbinding barrier and the binding free energy, and therefore favours guest unbinding both kinetically and thermodynamically. The analysis of different enthalpy contributions to the free energy shows that binding is enthalpically unfavourable and thus is an entropy-driven process, in agreement with the experimental data. Separating the enthalpy into the contributions from electrostatic, van der Waals, and bonded interactions in the force field shows that the unfavourable binding enthalpy is due to the bonded interactions being more favourable in the dissociated state, suggesting the presence of structural strain in the bound complex. Thus, the simulations provide microscopic explanations for the experimental findings and provide a possible route towards the targeted design of switchable nanocontainers with modified binding properties.

Olanzapine crystal symmetry originates in preformed centrosymmetric solute dimers

The symmetries of a crystal are notoriously uncorrelated to those of its constituent molecules. This symmetry breaking is typically thought to occur during crystallization. Here we demonstrate that one of the two symmetry elements of olanzapine crystals, an inversion centre, emerges in solute dimers extant in solution prior to crystallization. We combine time-resolved in situ scanning probe microscopy to monitor the crystal growth processes with all-atom molecular dynamics simulations. We show that crystals grow non-classically, predominantly by incorporation of centrosymmetric dimers. The growth rate of crystal layers exhibits a quadratic dependence on the solute concentration, characteristic of the second-order kinetics of the incorporation of dimers, which exist in equilibrium with a majority of monomers. We show that growth by dimers is preferred due to overwhelming accumulation of adsorbed dimers on the crystal surface, where it is complemented by dimerization and expedites dimer incorporation into growth sites.

Quantifying Protein-Protein Interactions in Molecular Simulations

Interactions among proteins, nucleic acids, and other macromolecules are essential for their biological functions and shape the physicochemcial properties of the crowded environments inside living cells. Binding interactions are commonly quantified by dissociation constants K_d, and both binding and nonbinding interactions are quantified by second osmotic virial coefficients B₂. As a measure of nonspecific binding and stickiness, B₂ is receiving renewed attention in the context of so-called liquid–liquid phase separation in protein and nucleic acid solutions. We show that K_d is fully determined by B₂ and the fraction of the dimer observed in molecular simulations of two proteins in a box. We derive two methods to calculate B₂. From molecular dynamics or Monte Carlo simulations using implicit solvents, we can determine B₂ from insertion and removal energies by applying Bennett’s acceptance ratio (BAR) method or the (binless) weighted histogram analysis method (WHAM). From simulations using implicit or explicit solvents, one can estimate B₂ from the probability that the two molecules are within a volume large enough to cover their range of interactions. We validate these methods for coarse-grained Monte Carlo simulations of three weakly binding proteins. Our estimates for K_d and B₂ allow us to separate out the contributions of nonbinding interactions to B₂. Comparison of calculated and measured values of K_d and B₂ can be used to (re-)parameterize and improve molecular force fields by calibrating specific affinities, overall stickiness, and nonbinding interactions. The accuracy and efficiency of K_d and B₂ calculations make them well suited for high-throughput studies of large interactomes.

Intrinsically disordered nuclear pore proteins show ideal-polymer morphologies and dynamics

In the nuclear pore complex, intrinsically disordered nuclear pore proteins (FG Nups) form a selective barrier for transport into and out of the cell nucleus, in a way that remains poorly understood. The collective FG Nup behavior has long been conceptualized either as a polymer brush, dominated by entropic and excluded-volume (repulsive) interactions, or as a hydrogel, dominated by cohesive (attractive) interactions between FG Nups. Here we compare mesoscale computational simulations with a wide range of experimental data to demonstrate that FG Nups are at the crossover point between these two regimes. Specifically, we find that repulsive and attractive interactions are balanced, resulting in morphologies and dynamics that are close to those of ideal polymer chains. We demonstrate that this property of FG Nups yields sufficient cohesion to seal the transport barrier, and yet maintains fast dynamics at the molecular scale, permitting the rapid polymer rearrangements needed for transport events.

Best Practices for Alchemical Free Energy Calculations [Article v1.0]

Alchemical free energy calculations are a useful tool for predicting free energy differences associated with the transfer of molecules from one environment to another. The hallmark of these methods is the use of "bridging" potential energy functions representing alchemical intermediate states that cannot exist as real chemical species. The data collected from these bridging alchemical thermodynamic states allows the efficient computation of transfer free energies (or differences in transfer free energies) with orders of magnitude less simulation time than simulating the transfer process directly. While these methods are highly flexible, care must be taken in avoiding common pitfalls to ensure that computed free energy differences can be robust and reproducible for the chosen force field, and that appropriate corrections are included to permit direct comparison with experimental data. In this paper, we review current best practices for several popular application domains of alchemical free energy calculations performed with equilibrium simulations, in particular relative and absolute small molecule binding free energy calculations to biomolecular targets.

Adsorption of Stearic Acid at the Iron Oxide/Oil Interface: Theory, Experiments, and Modeling

Improved friction performance is an important objective of equipment manufacturers for meeting improved energy efficiency demands. The addition of friction-reducing additives, or friction modifiers (FMs), to lubricants is a key part of the strategy. The performance of these additives is related to their surface activity and their ability to form adsorbed layers on the metal surface. However, the extent of surface coverage (mass per unit area) required for effective friction reduction is currently unknown. In this article, we show that full coverage is not necessary for significant friction reduction. We first highlight various features of surface adsorption that can influence the surface coverage, packing, and free energy of adsorption of organic FMs on iron oxide surfaces. Using stearic acid in heptane and hexadecane as model lubricant formulations, we employ a combination of experiments and molecular dynamics (MD) simulations to show how the dimerization of acid molecules in the bulk solvent and the crystallographic orientation of the surface modifies surface adsorption. In addition, we show that the solvent can strongly influence the adsorption kinetics, and MD simulations reveal that hexadecane tends to align on the surface, increasing the energy barrier for the adsorption of stearic acid to the surface. Furthermore, we present a combined approach using MD and molecular thermodynamic theory to calculate adsorption isotherms for stearic acid on iron oxide surfaces, which agrees well with experimental data obtained with a quartz crystal microbalance (QCM). Our results suggest that while the friction of systems lubricated with organic FMs decreases with increasing coverage, complete coverage of the surface is neither practically achievable nor necessary for effective friction reduction for the systems and conditions studied here.

Closely related, yet unique: Distinct homo- and heterodimerization patterns of G protein coupled chemokine receptors and their fine-tuning by cholesterol

Chemokine receptors, a subclass of G protein coupled receptors (GPCRs), play essential roles in the human immune system, they are involved in cancer metastasis as well as in HIV-infection. A plethora of studies show that homo- and heterodimers or even higher order oligomers of the chemokine receptors CXCR4, CCR5, and CCR2 modulate receptor function. In addition, membrane cholesterol affects chemokine receptor activity. However, structural information about homo- and heterodimers formed by chemokine receptors and their interplay with cholesterol is limited. Here, we report homo- and heterodimer configurations of the chemokine receptors CXCR4, CCR5, and CCR2 at atomistic detail, as obtained from thousands of molecular dynamics simulations. The observed homodimerization patterns were similar for the closely related CC chemokine receptors, yet they differed significantly between the CC receptors and CXCR4. Despite their high sequence identity, cholesterol modulated the CC homodimer interfaces in a subtype-specific manner. Chemokine receptor heterodimers display distinct dimerization patterns for CXCR4/CCR5 and CXCR4/CCR2. Furthermore, associations between CXCR4 and CCR5 reveal an increased cholesterol-sensitivity as compared to CXCR4/CCR2 heterodimerization patterns. This work provides a first comprehensive structural overview over the complex interaction network between chemokine receptors and indicates how heterodimerization and the interaction with the membrane environment diversifies the function of closely related GPCRs.

Classical Molecular Dynamics with Mobile Protons

An important limitation of standard classical molecular dynamics simulations is the inability to make or break chemical bonds. This restricts severely our ability to study processes that involve even the simplest of chemical reactions, the transfer of a proton. Existing approaches for allowing proton transfer in the context of classical mechanics are rather cumbersome and have not achieved widespread use and routine status. Here we reconsider the combination of molecular dynamics with periodic stochastic proton hops. To ensure computational efficiency, we propose a non-Boltzmann acceptance criterion that is heuristically adjusted to maintain the correct or desirable thermodynamic equilibria between different protonation states and proton transfer rates. Parameters are proposed for hydronium, Asp, Glu, and His. The algorithm is implemented in the program CHARMM and tested on proton diffusion in bulk water and carbon nanotubes and on proton conductance in the gramicidin A channel. Using hopping parameters determined from proton diffusion in bulk water, the model reproduces the enhanced proton diffusivity in carbon nanotubes and gives a reasonable estimate of the proton conductance in gramicidin A.

Self-association of a highly charged arginine-rich cell-penetrating peptide

Small-angle X-ray scattering (SAXS) measurements reveal a striking difference in intermolecular interactions between two short highly charged peptides-deca-arginine (R10) and deca-lysine (K10). Comparison of SAXS curves at high and low salt concentration shows that R10 self-associates, while interactions between K10 chains are purely repulsive. The self-association of R10 is stronger at lower ionic strengths, indicating that the attraction between R10 molecules has an important electrostatic component. SAXS data are complemented by NMR measurements and potentials of mean force between the peptides, calculated by means of umbrella-sampling molecular dynamics (MD) simulations. All-atom MD simulations elucidate the origin of the R10-R10 attraction by providing structural information on the dimeric state. The last two C-terminal residues of R10 constitute an adhesive patch formed by stacking of the side chains of two arginine residues and by salt bridges formed between the like-charge ion pair and the C-terminal carboxyl groups. A statistical analysis of the Protein Data Bank reveals that this mode of interaction is a common feature in proteins.

Structural and Spectroscopic Properties of Assemblies of Self-Replicating Peptide Macrocycles

Self-replication at the molecular level is often seen as essential to the early origins of life. Recently a mechanism of self-replication has been discovered in which replicator self-assembly drives the process. We have studied one of the examples of such self-assembling self-replicating molecules to a high level of structural detail using a combination of computational and spectroscopic techniques. Molecular Dynamics simulations of self-assembled stacks of peptide-derived replicators provide insights into the structural characteristics of the system and serve as the basis for semiempirical calculations of the UV-vis, circular dichroism (CD) and infrared (IR) absorption spectra that reflect the chiral organization and peptide secondary structure of the stacks. Two proposed structural models are tested by comparing calculated spectra to experimental data from electron microscopy, CD and IR spectroscopy, resulting in a better insight into the specific supramolecular interactions that lead to self-replication. Specifically, we find a cooperative self-assembly process in which β-sheet formation leads to well-organized structures, while also the aromatic core of the macrocycles plays an important role in the stability of the resulting fibers.

Atomistic fingerprint of hyaluronan-CD44 binding

Hyaluronan is a polyanionic, megadalton-scale polysaccharide, which initiates cell signaling by interacting with several receptor proteins including CD44 involved in cell-cell interactions and cell adhesion. Previous studies of the CD44 hyaluronan binding domain have identified multiple widespread residues to be responsible for its recognition capacity. In contrast, the X-ray structural characterization of CD44 has revealed a single binding mode associated with interactions that involve just a fraction of these residues. In this study, we show through atomistic molecular dynamics simulations that hyaluronan can bind CD44 with three topographically different binding modes that in unison define an interaction fingerprint, thus providing a plausible explanation for the disagreement between the earlier studies. Our results confirm that the known crystallographic mode is the strongest of the three binding modes. The other two modes represent metastable configurations that are readily available in the initial stages of the binding, and they are also the most frequently observed modes in our unbiased simulations. We further discuss how CD44, fostered by the weaker binding modes, diffuses along HA when attached. This 1D diffusion combined with the constrained relative orientation of the diffusing proteins is likely to influence the aggregation kinetics of CD44. Importantly, CD44 aggregation has been suggested to be a possible mechanism in CD44-mediated signaling.

Hyaluronan is a natural sugar polymer in our bodies. Besides acting as a space-filling agent for example in multiple connective tissues, it also functions as a cellular cue in cancer and inflammation. Our tissues sense hyaluronan through receptors—proteins that sit at the surface of cells and grab the molecules they are expected to recognize. Although the knowledge associated with hyaluronan and its receptors is constantly accumulating, the molecular-level insight is largely missing or incomplete due to the lack of techniques able to probe the dynamics of protein–carbohydrate interactions with sufficiently high resolution. In this work, we characterize the binding of hyaluronan to its receptor CD44 with atomistic precision. We achieve this level of precision by employing atomistic molecular dynamics simulations. This computational technique allows one to follow the movement of atoms of a virtual system at scales beyond the resolution of any experimental technique. Our work specifically focuses on the different stages of hyaluronan–CD44 binding, and we observe the process to involve three different binding modes, making it more versatile than previously thought. Our insights, therefore, promote the understanding of the interplay between hyaluronan and HA, thereby fostering development of new drugs or inhibitors to malignancies, such as cancer metastasis.

Quantitative Characterization of the Binding and Unbinding of Millimolar Drug Fragments with Molecular Dynamics Simulations

A quantitative characterization of the binding properties of drug fragments to a target protein is an important component of a fragment-based drug discovery program. Fragments typically have a weak binding affinity, however, making it challenging to experimentally characterize key binding properties, including binding sites, poses, and affinities. Direct simulation of the binding equilibrium by molecular dynamics (MD) simulations can provide a computational route to characterize fragment binding, but this approach is so computationally intensive that it has thus far remained relatively unexplored. Here, we perform MD simulations of sufficient length to observe several different fragments spontaneously and repeatedly bind to and unbind from the protein FKBP, allowing the binding affinities, on- and off-rates, and relative occupancies of alternative binding sites and alternative poses within each binding site to be estimated, thereby illustrating the potential of long time scale MD as a quantitative tool for fragment-based drug discovery. The data from the long time scale fragment binding simulations reported here also provide a useful benchmark for testing alternative computational methods aimed at characterizing fragment binding properties. As an example, we calculated binding affinities for the same fragments using a standard free energy perturbation approach and found that the values agreed with those obtained from the fragment binding simulations within statistical error.

Effect of Sodium and Chloride Binding on a Lecithin Bilayer. A Molecular Dynamics Study

The effect of ion binding on the structural, mechanical, dynamic and electrostatic properties of a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer in a 0.5 M aqueous NaCl solution is investigated using classical atomistic molecular dynamics simulation with different force-field descriptions for ion-ion and ion-lipid interactions. Most importantly, the repulsive Lennard-Jones parameters for the latter were modified, such that approximately similar binding of cations and anions to the lipid membrane is achieved. This was done to qualitatively improve the apparent ion-lipid binding constants obtained from simulations with the original force field (Berger lipids and GROMOS87 ions in combination with the SPC water model) in comparison to experimental data. Furthermore, various parameters characterizing membrane structure, elasticity, order and dynamics are analyzed. It is found that ion binding as observed in simulations involving the modified in comparison to the original force-field description leads to: (i) a smaller salt-induced change in the area per lipid, which is in closer agreement with the experiment; (ii) a decrease in the area compressibility and bilayer thickness to values comparable to a bilayer in pure water; (iii) lipid deuterium order parameters and lipid diffusion coefficients on nanosecond timescales that are very similar to the values for a membrane in pure water. In general, salt effects on the structural properties of a POPC bilayer in an aqueous sodium-chloride solution appear to be reproduced reasonably well by the new force-field description. An analysis of membrane-membrane disjoining pressure suggests that the smaller salt-induced change in area per lipid induced by the new force-field description is not due to the alteration of membrane-associated net charge, but must rather be understood as a consequence of ion-specific effects on the arrangement of lipid molecules.

Free-energy patterns in inclusion complexes: the relevance of non-included moieties in the stability constants

A MD/PMF-based procedure is designed for quantification of the interaction and respective components, guiding complex formation in water between β-CD and several naphthalene derivatives, highlighting the relevance of substituents.

Convergence and Sampling in Determining Free Energy Landscapes for Membrane Protein Association

Potential of mean force (PMF) calculations are used to characterize the free energy landscape of protein–lipid and protein–protein association within membranes. Coarse-grained simulations allow binding free energies to be determined with reasonable statistical error. This accuracy relies on defining a good collective variable to describe the binding and unbinding transitions, and upon criteria for assessing the convergence of the simulation toward representative equilibrium sampling. As examples, we calculate protein–lipid binding PMFs for ANT/cardiolipin and Kir2.2/PIP₂, using umbrella sampling on a distance coordinate. These highlight the importance of replica exchange between windows for convergence. The use of two independent sets of simulations, initiated from bound and unbound states, provide strong evidence for simulation convergence. For a model protein–protein interaction within a membrane, center-of-mass distance is shown to be a poor collective variable for describing transmembrane helix–helix dimerization. Instead, we employ an alternative intermolecular distance matrix RMS (D_RMS) coordinate to obtain converged PMFs for the association of the glycophorin transmembrane domain. While the coarse-grained force field gives a reasonable K_d for dimerization, the majority of the bound population is revealed to be in a near-native conformation. Thus, the combination of a refined reaction coordinate with improved sampling reveals previously unnoticed complexities of the dimerization free energy landscape. We propose the use of replica-exchange umbrella sampling starting from different initial conditions as a robust approach for calculation of the binding energies in membrane simulations.

Dynamic Cholesterol-Conditioned Dimerization of the G Protein Coupled Chemokine Receptor Type 4

G protein coupled receptors (GPCRs) allow for the transmission of signals across biological membranes. For a number of GPCRs, this signaling was shown to be coupled to prior dimerization of the receptor. The chemokine receptor type 4 (CXCR4) was reported before to form dimers and their functionality was shown to depend on membrane cholesterol. Here, we address the dimerization pattern of CXCR4 in pure phospholipid bilayers and in cholesterol-rich membranes. Using ensembles of molecular dynamics simulations, we show that CXCR4 dimerizes promiscuously in phospholipid membranes. Addition of cholesterol dramatically affects the dimerization pattern: cholesterol binding largely abolishes the preferred dimer motif observed for pure phospholipid bilayers formed mainly by transmembrane helices 1 and 7 (TM1/TM5-7) at the dimer interface. In turn, the symmetric TM3,4/TM3,4 interface is enabled first by intercalating cholesterol molecules. These data provide a molecular basis for the modulation of GPCR activity by its lipid environment.

An Amphotericin B Derivative Equally Potent to Amphotericin B and with Increased Safety

Amphotericin B is the most potent antimycotic known to date. However due to its large collateral toxicity, its use, although long standing, had been limited. Many attempts have been made to produce derivatives with reduced collateral damage. The molecular mechanism of polyene has also been closely studied for this purpose and understanding it would contribute to the development of safe derivatives. Our study examined polyene action, including chemical synthesis, electrophysiology, pharmacology, toxicology and molecular dynamics. The results were used to support a novel Amphotericin B derivative with increased selectivity: L-histidine methyl ester of Amphotericin B. We found that this derivative has the same form of action as Amphotericin B, i.e. pore formation in the cell membrane. Its reduced dimerization in solution, when compared to Amphotericin B, is at least partially responsible for its increased selectivity. Here we also present the results of preclinical tests, which show that the derivative is just as potent as Amphotericin B and has increased safety.

A nucleotide-level coarse-grained model of RNA

We present a new, nucleotide-level model for RNA, oxRNA, based on the coarse-graining methodology recently developed for the oxDNA model of DNA. The model is designed to reproduce structural, mechanical, and thermodynamic properties of RNA, and the coarse-graining level aims to retain the relevant physics for RNA hybridization and the structure of single- and double-stranded RNA. In order to explore its strengths and weaknesses, we test the model in a range of nanotechnological and biological settings. Applications explored include the folding thermodynamics of a pseudoknot, the formation of a kissing loop complex, the structure of a hexagonal RNA nanoring, and the unzipping of a hairpin motif. We argue that the model can be used for efficient simulations of the structure of systems with thousands of base pairs, and for the assembly of systems of up to hundreds of base pairs. The source code implementing the model is released for public use.

Using interpolation for fast and accurate calculation of ion-ion interactions

We perform extensive molecular dynamics (MD) simulations between pairs of ions of various diameters (2-5.5 Å in increments of 0.5 Å) and charge (+1 or -1) interacting in explicit water (TIP3P) under ambient conditions. We extract their potentials of mean force (PMFs). We develop an interpolation scheme, called i-PMF, that is capable of capturing the full set of PMFs for arbitrary combinations of ion sizes ranging from 2 to 5.5 Å. The advantage of the interpolation process is computational cost. Whereas it can take 100 h to simulate each PMF by MD, we can compute an equivalently accurate i-PMF in seconds. This process may be useful for rapid and accurate calculation of the strengths of salt bridges and the effects of bridging waters in biomolecular simulations. We also find that our data is consistent with Collins' "law of matching affinities" of ion solubilities: small-small or large-large ion pairs are poorly soluble in water, whereas small-large are highly soluble.

Coarse-graining DNA for simulations of DNA nanotechnology

To simulate long time and length scale processes involving DNA it is necessary to use a coarse-grained description. Here we provide an overview of different approaches to such coarse-graining, focussing on those at the nucleotide level that allow the self-assembly processes associated with DNA nanotechnology to be studied. OxDNA, our recently-developed coarse-grained DNA model, is particularly suited to this task, and has opened up this field to systematic study by simulations. We illustrate some of the range of DNA nanotechnology systems to which the model is being applied, as well as the insights it can provide into fundamental biophysical properties of DNA.