Exploring Multidimensional Free Energy Landscapes Using Time-Dependent Biases on Collective Variables

Interaction between SARS-CoV-2 spike glycoprotein and human skin models: a molecular dynamics study

The possibility of contamination of human skin by infectious virions plays an important role in indirect transmission of respiratory viruses but little is known about the fundamental physico-chemical aspects of the virus-skin interactions. In the case of coronaviruses, the interaction with surfaces (including the skin surface) is mediated by their large glycoprotein spikes that protrude from (and cover) the viral envelope. Here, we perform all atomic simulations between the SARS-CoV-2 spike glycoprotein and human skin models. We consider an "oily" skin covered by sebum and a "clean" skin exposing the stratum corneum. The simulations show that the spike tries to maximize the contacts with stratum corneum lipids, particularly ceramides, with substantial hydrogen bonding. In the case of "oily" skin, the spike is able to retain its structure, orientation and hydration over sebum with little interaction with sebum components. Comparison of these results with our previous simulations of the interaction of SARS-CoV-2 spike with hydrophilic and hydrophobic solid surfaces, suggests that the "soft" or "hard" nature of the surface plays an essential role in the interaction of the spike protein with materials.

Toward Convergence in Free Energy Calculations for Protein Conformational Changes: A Case Study on the Thin Gate of Mhp1 Transporter

It has been challenging to obtain reliable free energies for protein conformational changes from all-atom molecular dynamics simulations, despite the availability of many enhanced sampling techniques. To alleviate the difficulties associated with the enormous complexity of the conformational space, here we propose a few practical strategies for such calculations, including (1) a stringent method to examine convergence by comparing independent simulations starting from different initial coordinates, (2) adoption of multistep schemes in which the complete conformational change consists of multiple transition steps, each sampled using a distinct reaction coordinate, and (3) application of boundary restraints to simplify the conformational space. We demonstrate these strategies on the conformational changes between the outward-facing and outward-occluded states of the Mhp1 membrane transporter, obtaining the equilibrium thermodynamics of the relevant metastable states, the kinetic rates between these states, and the reactive trajectories that reveal the atomic details of spontaneous transitions. Our approaches thus promise convergent and reliable calculations to examine intuition-based hypotheses and to eventually elucidate the underlying molecular mechanisms of reversible conformational changes in complex protein systems.

Molecular dynamics simulations of doxorubicin in sphingomyelin-based lipid membranes

Doxorubicin (DOX) is one of the most efficient antitumor drugs employed in numerous cancer therapies. Its incorporation into lipid-based nanocarriers, such as liposomes, improves the drug targeting into tumor cells and reduces drug side effects. The carriers' lipid composition is expected to affect the interactions of DOX and its partitioning into liposomal membranes. To get a rational insight into this aspect and determine promising lipid compositions, we use numerical simulations, which provide unique information on DOX-membrane interactions at the atomic level of resolution. In particular, we combine classical molecular dynamics simulations and free energy calculations to elucidate the mechanism of penetration of a protonated Doxorubicin molecule (DOX⁺) into potential liposome membranes, here modeled as lipid bilayers based on mixtures of phosphatidylcholine (PC), sphingomyelin (SM) and cholesterol lipid molecules, of different compositions and lipid phases. Moreover, we analyze DOX⁺ partitioning into relevant regions of SM-based lipid bilayer systems using a combination of free energy methods. Our results show that DOX⁺ penetration and partitioning are facilitated into less tightly packed SM-based membranes and are dependent on lipid composition. This work paves the way to further investigations of optimal formulations for lipid-based carriers, such as those associated with pH-responsive membranes.

Correlated Response of Protein Side-Chain Fluctuations and Conformational Entropy to Ligand Binding

The heterogeneous fast side-chain dynamics of proteins plays crucial roles in molecular recognition and binding. Site-specific NMR experiments quantify these motions by measuring the model-free order parameter (O_axis²) on a scale of 0 (most flexible) to 1 (least flexible) for each methyl-containing residue of proteins. Here, we have examined ligand-induced variations in the fast side-chain dynamics and conformational entropy of calmodulin (CaM) using five different CaM-peptide complexes. O_axis² of CaM in the ligand-free (O_axis,U²) and ligand-bound (O_axis,B²) states are calculated from molecular dynamics trajectories and conformational energy surfaces obtained using the adaptive biasing force (ABF) method. ΔO_axis² = O_axis,B² - O_axis,U² follows a Gaussian-like unimodal distribution whose second moment is a potential indicator of the binding affinity of these complexes. The probability for the binding-induced O_axis,U² → O_axis,B² transition decreases with increasing magnitude of ΔO_axis², indicating that large flexibility changes are improbable for side chains of CaM after ligand binding. A linear correlation established between ΔO_axis² and the conformational entropy change of the protein makes possible the determination of the conformational entropy of binding of protein-ligand complexes. The results not only underscore the functional importance of fast side-chain fluctuations but also highlight key motional and thermodynamic correlates of protein-ligand binding.

Folding and Intrinsic Disorder of the Receptor Tyrosine Kinase KIT Insert Domain Seen by Conventional Molecular Dynamics Simulations

The kinase insert domain (KID) of RTK KIT is the key recruitment region for downstream signalling proteins. KID, studied by molecular dynamics simulations as a cleaved polypeptide and as a native domain fused to KIT, showed intrinsic disorder represented by a set of heterogeneous conformations. The accurate atomistic models showed that the helical fold of KID is mainly sequence dependent. However, the reduced fold of the native KID suggests that its folding is allosterically controlled by the kinase domain. The tertiary structure of KID represents a compact array of highly variable α- and 3₁₀-helices linked by flexible loops playing a principal role in the conformational diversity. The helically folded KID retains a collapsed globule-like shape due to non-covalent interactions associated in a ternary hydrophobic core. The free energy landscapes constructed from first principles-the size, the measure of the average distance between the conformations, the amount of helices and the solvent-accessible surface area-describe the KID disorder through a collection of minima (wells), providing a direct evaluation of conformational ensembles. We found that the cleaved KID simulated with restricted N- and C-ends better reproduces the native KID than the isolated polypeptide. We suggest that a cyclic, generic KID would be best suited for future studies of KID f post-transduction effects.

Optimizing String Method's Reproducibility Using Generalized Solute Tempering Replica Exchange

Obtaining accurate and reproducible free energies from molecular simulations is somewhat tricky due to incomplete knowledge of crucial slow degrees of freedom leading to hidden barriers that can stymie sampling. Employing a sufficiently large number of collective variables (CV) and ensuring ergodic sampling in orthogonal CV space, perhaps via tempering methods, can reduce these issues to some extent. For complex systems with high-dimensional free energy landscapes, both these approaches become computationally expensive. For high-dimensional landscapes, efficient exploration can be enabled by using temperature-accelerated MD (TAMD) and identification and characterization of minimum free energy pathways connecting minima can be found by using the string method (SM). Both TAMD and SM use mean-force estimates from finite MD simulations and are thus susceptible to sampling restrictions from hidden variables. A recent development in parallel tempering methods, "generalized replica exchange solute tempering" (gREST), can enhance sampling at a reasonable computational cost with its flexibility to target very specific "solutes" which can include arbitrary independent variables. Considering the advantages of both methods, we implement gREST-enabled TAMD and SM. By considering two different collective variable representations of the pentapeptide neurotransmitter met-enkephalin, we show that both gREST-enabled TAMD and SM yield more accurate and reproducible free energy predictions than TAMD and SM alone. Given the moderate computational cost of gREST compared with other replica-exchange methods, gREST-enabled SM represents a more attractive method for characterizing free energy minima and pathways among them for a large variety of systems.

Free energy landscape of RNA binding dynamics in start codon recognition by eukaryotic ribosomal pre-initiation complex

Specific interaction between the start codon, 5'-AUG-3', and the anticodon, 5'-CAU-3', ensures accurate initiation of translation. Recent studies show that several near-cognate start codons (e.g. GUG and CUG) can play a role in initiating translation in eukaryotes. However, the mechanism allowing initiation through mismatched base-pairs at the ribosomal decoding site is still unclear at an atomic level. In this work, we propose an extended simulation-based method to evaluate free energy profiles, through computing the distance between each base-pair of the triplet interactions involved in recognition of start codons in eukaryotic translation pre-initiation complex. Our method provides not only the free energy penalty for mismatched start codons relative to the AUG start codon, but also the preferred pathways of transitions between bound and unbound states, which has not been described by previous studies. To verify the method, the binding dynamics of cognate (AUG) and near-cognate start codons (CUG and GUG) were simulated. Evaluated free energy profiles agree with experimentally observed changes in initiation frequencies from respective codons. This work proposes for the first time how a G:U mismatch at the first position of codon (GUG)-anticodon base-pairs destabilizes the accommodation in the initiating eukaryotic ribosome and how initiation at a CUG codon is nearly as strong as, or sometimes stronger than, that at a GUG codon. Our method is expected to be applied to study the affinity changes for various mismatched base-pairs.

Effects of Cholesterol on the Partitioning of a Drug Molecule in Lipid Bilayers

Drug molecules either bind to membrane-bound targets or permeate through cell membranes to reach intracellular targets, and hence, their membrane partition and permeation are of great importance. Here, we studied the effects of cholesterol on the partition of amantadine, an antiflu drug molecule, into 1,2-dimyristoyl-sn-glycero-3-phosphocholine bilayers using molecular dynamics simulations. The membrane partition of amantadine is sensitive to the cholesterol mole fraction (x_chol). In the absence of cholesterol, amantadine is stably bound in membranes, but at x_chol = 32%, it can escape to the aqueous phase, in agreement with recent experiments. The reduced membrane partition of amantadine at a high cholesterol content is mainly due to the perturbation of the bilayer structure and dynamics. Surrounding lipids stabilize amantadine by having their tails wrapped around the drug molecule, and this ability is compromised when cholesterol is present to increase the order in lipid tails. The atomic details on interactions with lipids and perturbations by cholesterol revealed here provide insight into membrane partition and delivery of drug molecules to their targets.

Steric Hindrance in Metal Coordination Drives the Separation of Pyridine Regioisomers Using Rhodium(II)-Based Metal-Organic Polyhedra

The physicochemical similarity of isomers makes their chemical separation through conventional techniques energy intensive. Herein, we report that, instead of using traditional encapsulation-driven processes, steric hindrance in metal coordination on the outer surface of Rh^II -based metal-organic polyhedra (Rh-MOPs) can be used to separate pyridine-based regioisomers via liquid-liquid extraction. Through molecular dynamics simulations and wet experiments, we discovered that the capacity of pyridines to coordinatively bind to Rh-MOPs is determined by the positions of the pyridine substituents relative to the pyridine nitrogen and is influenced by steric hindrance. Thus, we exploited the differential solubility of bound and non-bound pyridine regioisomers to engineer liquid-liquid self-sorting systems. As a proof of concept, we separated four different equimolecular mixtures of regioisomers, including a mixture of the industrially relevant compounds 2-chloropyridine and 3-chloropyridine, isolating highly pure compounds in all cases.

An Asymmetric Mechanism in a Symmetric Molecular Machine

The design of molecular architectures exhibiting functional motions is a promising area for disruptive technological development. Toward this goal, rotaxanes and catenanes, which undergo relative motions of their subunits in response to external stimuli, are prime candidates. Here, we report on the computational analysis of the contraction/extension of a bistable [c2]daisy chain rotaxane. Using free-energy calculations and transition path optimizations, we explore the free-energy landscape governing the functional motions of a prototypical molecular machine with atomic resolution. The calculations reveal a sequential mechanism in which the asynchronous gliding of each ring is preferred over the concerted movement. Analysis of the underlying free-energy surface indicates that the formation of partially rearranged intermediates entails crossing of much smaller barriers. Our findings illustrate an important design principle for molecular machines, namely that efficient exploitation of thermal fluctuations may be realized by breaking down the large-scale functional motions into smaller steps.

Microscopic Interactions of Melatonin, Serotonin and Tryptophan with Zwitterionic Phospholipid Membranes

The interactions at the atomic level between small molecules and the main components of cellular plasma membranes are crucial for elucidating the mechanisms allowing for the entrance of such small species inside the cell. We have performed molecular dynamics and metadynamics simulations of tryptophan, serotonin, and melatonin at the interface of zwitterionic phospholipid bilayers. In this work, we will review recent computer simulation developments and report microscopic properties, such as the area per lipid and thickness of the membranes, atomic radial distribution functions, angular orientations, and free energy landscapes of small molecule binding to the membrane. Cholesterol affects the behaviour of the small molecules, which are mainly buried in the interfacial regions. We have observed a competition between the binding of small molecules to phospholipids and cholesterol through lipidic hydrogen-bonds. Free energy barriers that are associated to translational and orientational changes of melatonin have been found to be between 10-20 kJ/mol for distances of 1 nm between melatonin and the center of the membrane. Corresponding barriers for tryptophan and serotonin that are obtained from reversible work methods are of the order of 10 kJ/mol and reveal strong hydrogen bonding between such species and specific phospholipid sites. The diffusion of tryptophan and melatonin is of the order of 10-7 cm2/s for the cholesterol-free and cholesterol-rich setups.

Benchmarking Force Field and the ANI Neural Network Potentials for the Torsional Potential Energy Surface of Biaryl Drug Fragments

Many drug molecules contain biaryl fragments, resulting in a torsional barrier corresponding to rotation around the bond linking the aryls. The potential energy surfaces of these torsions vary significantly because of steric and electronic effects, ultimately affecting the relative stability of the molecular conformations in the protein-bound and solution states. Simulations of protein-ligand binding require accurate computational models to represent the intramolecular interactions to provide accurate predictions of the structure and dynamics of binding. In this article, we compare four force fields [generalized AMBER force field (GAFF), open force field (OpenFF), CHARMM general force field (CGenFF), optimized potentials for liquid simulations (OPLS)] and two neural network potentials (ANI-2x and ANI-1ccx) for their ability to predict the torsional potential energy surfaces of 88 biaryls extracted from drug fragments. The root mean square deviation (rmsd) over the full potential energy surface and the mean absolute deviation of the torsion rotational barrier height (MADB) relative to high-level ab initio reference data (CCSD(T1)*) were used as the measure of accuracy. Uncertainties in these metrics due to the composition of the data set were estimated using bootstrap analysis. In comparison to high-level ab initio data, ANI-1ccx was most accurate for predicting the barrier height (rmsd: 0.5 ± 0.0 kcal/mol, MADB: 0.8 ± 0.1 kcal/mol), followed closely by ANI-2x (rmsd: 0.5 ± 0.0 kcal/mol, MADB: 1.0 ± 0.2 kcal/mol), then CGenFF (rmsd: 0.8 ± 0.1 kcal/mol, MADB: 1.3 ± 0.1 kcal/mol) and OpenFF (rmsd: 0.7 ± 0.1 kcal/mol, MADB: 1.3 ± 0.1 kcal/mol), then GAFF (rmsd: 1.2 ± 0.2 kcal/mol, MADB: 2.6 ± 0.5 kcal/mol), and finally OPLS (rmsd: 3.6 ± 0.3 kcal/mol, MADB: 3.6 ± 0.3 kcal/mol). Significantly, the neural network potentials (NNPs) are systematically more accurate and more reliable than any of the force fields. As a practical example, the NNP/molecular mechanics method was used to simulate the isomerization of ozanimod, a drug used for multiple sclerosis. Multinanosecond molecular dynamics (MD) simulations in an explicit aqueous solvent were performed, as well as umbrella sampling and adaptive biasing force-enhanced sampling techniques. The rate constant for this isomerization calculated using transition state theory was 4.30 × 10^-1 ns^-1, which is consistent with direct MD simulations.

Long-lasting Salt Bridges Provide the Anchoring Mechanism of Oncogenic Kirsten Rat Sarcoma Proteins at Cell Membranes

RAS proteins work as GDP-GTP binary switches and regulate cytoplasmic signaling networks that are able to control several cellular processes, playing an essential role in signal transduction pathways involved in cell growth, differentiation, and survival, so that overacting RAS signaling can lead to cancer. One of the hardest challenges to face is the design of mutation-selective therapeutic strategies. In this work, a G12D-mutated farnesylated GTP-bound Kirsten RAt sarcoma (KRAS) protein has been simulated at the interface of a DOPC/DOPS/cholesterol model anionic cell membrane. A specific long-lasting salt bridge connection between farnesyl and the hypervariable region of the protein has been identified as the main mechanism responsible for the binding of oncogenic farnesylated KRAS-4B to the cell membrane. Free-energy landscapes allowed us to characterize local and global minima of KRAS-4B binding to the cell membrane, revealing the main pathways between anchored and released states.

Thermodynamics and Mechanism of the Membrane Permeation of Hv1 Channel Blockers

The voltage-gated proton channel Hv1 mediates efflux of protons from the cell. Hv1 integrally contributes to various physiological processes including pH homeostasis and the respiratory burst of phagocytes. Inhibition of Hv1 may provide therapeutic avenues for the treatment of inflammatory diseases, breast cancer, and ischemic brain damage. In this work, we investigate two prototypical Hv1 inhibitors, 2-guanidinobenzimidazole (2GBI), and 5-chloro-2-guanidinobenzimidazole (GBIC), from an experimentally screened class of guanidine derivatives. Both compounds block proton conduction by binding the same site located on the intracellular side of the channel. However, when added to the extracellular medium, the compounds strongly differ in their ability to inhibit proton conduction, suggesting substantial differences in membrane permeability. Here, we compute the potential of mean force for each compound to permeate through the membrane using atomistic molecular dynamics simulations with the adaptive biasing force method. Our results rationalize the putative distinction between these two blockers with respect to their abilities to permeate the cellular membrane.

Interaction Mechanism of the Germination Stimulants Karrikins and Their Receptor ShKAI2iB

The significance of karrikins (KARs) in plant physiology opens a door for their application in the agricultural production. As the first event of the whole signaling pathway, the binding of smoke-derived signal molecules KARs to the receptor protein KAI2 triggers the germination of the primary dormant seeds of all angiosperms, not just the "fire-prone" taxa. In the present study, all-atom molecular dynamics simulations, along with the accurate estimation of the ligand-receptor binding free energy, were used to investigate the atomic level interaction of all the members of the KARs family (from KAR₁ to KAR₆) with the receptor ShKAI2iB, an intermediate-evolving KAI2 from Striga hermonthica. The calculated binding energy value of KAR₁ to ShKAI2iB, -5.64 kcal/mol, is in good agreement with the available experimental data, -5.67 kcal/mol. The further analysis of the detailed interaction between each KAR and the protein reveals the primary reasons for the difference of the affinity of the diverse ligands with the receptor and displays the regional characteristics of the protein's active site. Our research will not only provide clues for the study of equivalent endogenous phytohormone, but also contribute to the development of synthetic germinating chemicals.

Gamma estimator of Jarzynski equality for recovering binding energies from noisy dynamic data sets

A fundamental problem in thermodynamics is the recovery of macroscopic equilibrated interaction energies from experimentally measured single-molecular interactions. The Jarzynski equality forms a theoretical basis in recovering the free energy difference between two states from exponentially averaged work performed to switch the states. In practice, the exponentially averaged work value is estimated as the mean of finite samples. Numerical simulations have shown that samples having thousands of measurements are not large enough for the mean to converge when the fluctuation of external work is above 4 k_BT, which is easily observable in biomolecular interactions. We report the first example of a statistical gamma work distribution applied to single molecule pulling experiments. The Gibbs free energy of surface adsorption can be accurately evaluated even for a small sample size. The values obtained are comparable to those derived from multi-parametric surface plasmon resonance measurements and molecular dynamics simulations.

Measuring interaction energies from experimentally measured single-molecular interactions is challenging. Here, the authors report a gamma work distribution applied to single molecule pulling events for estimating peptide absorption free energy.

Computer simulations of the interaction between SARS-CoV-2 spike glycoprotein and different surfaces

A prominent feature of coronaviruses is the presence of a large glycoprotein spike protruding from a lipidic membrane. This glycoprotein spike determines the interaction of coronaviruses with the environment and the host. In this paper, we perform all atomic molecular dynamics simulations of the interaction between the SARS-CoV-2 trimeric glycoprotein spike and surfaces of materials. We considered a material with high hydrogen bonding capacity (cellulose) and a material capable of strong hydrophobic interactions (graphite). Initially, the spike adsorbs to both surfaces through essentially the same residues belonging to the receptor binding subunit of its three monomers. Adsorption onto cellulose stabilizes in this configuration, with the help of a large number of hydrogen bonds developed between cellulose and the three receptor-binding domains of the glycoprotein spike. In the case of adsorption onto graphite, the initial adsorption configuration is not stable and the surface induces a substantial deformation of the glycoprotein spike with a large number of adsorbed residues not pertaining to the binding subunits of the spike monomers.

A novel multiscale scheme to accelerate atomistic simulations of bio-macromolecules by adaptively driving coarse-grained coordinates

All-atom molecular dynamics (MD) simulations of bio-macromolecules can yield relatively accurate results while suffering from the limitation of insufficient conformational sampling. On the other hand, the coarse-grained (CG) MD simulations efficiently accelerate conformational changes in biomolecules but lose atomistic details and accuracy. Here, we propose a novel multiscale simulation method called the adaptively driving multiscale simulation (ADMS)-it efficiently accelerates biomolecular dynamics by adaptively driving virtual CG atoms on the fly while maintaining the atomistic details and focusing on important conformations of the original system with irrelevant conformations rarely sampled. Herein, the "adaptive driving" is based on the short-time-averaging response of the system (i.e., an approximate free energy surface of the original system), without requiring the construction of the CG force field. We apply the ADMS to two peptides (deca-alanine and Ace-GGPGGG-Nme) and one small protein (HP35) as illustrations. The simulations show that the ADMS not only efficiently captures important conformational states of biomolecules and drives fast interstate transitions but also yields, although it might be in part, reliable protein folding pathways. Remarkably, a ∼100-ns explicit-solvent ADMS trajectory of HP35 with three CG atoms realizes folding and unfolding repeatedly and captures the important states comparable to those from a 398-µs standard all-atom MD simulation.

Structure and Function of the T4 Spackle Protein Gp61.3

The bacteriophage T4 genome contains two genes that code for proteins with lysozyme activity—e and 5. Gene e encodes the well-known T4 lysozyme (commonly called T4L) that functions to break the peptidoglycan layer late in the infection cycle, which is required for liberating newly assembled phage progeny. Gene product 5 (gp5) is the tail-associated lysozyme, a component of the phage particle. It forms a spike at the tip of the tail tube and functions to pierce the outer membrane of the Escherichia coli host cell after the phage has attached to the cell surface. Gp5 contains a T4L-like lysozyme domain that locally digests the peptidoglycan layer upon infection. The T4 Spackle protein (encoded by gene 61.3) has been thought to play a role in the inhibition of gp5 lysozyme activity and, as a consequence, in making cells infected by bacteriophage T4 resistant to later infection by T4 and closely related phages. Here we show that (1) gp61.3 is secreted into the periplasm where its N-terminal periplasm-targeting peptide is cleaved off; (2) gp61.3 forms a 1:1 complex with the lysozyme domain of gp5 (gp5Lys); (3) gp61.3 selectively inhibits the activity of gp5, but not that of T4L; (4) overexpression of gp5 causes cell lysis. We also report a crystal structure of the gp61.3-gp5Lys complex that demonstrates that unlike other known lysozyme inhibitors, gp61.3 does not interact with the active site cleft. Instead, it forms a “wall” that blocks access of an extended polysaccharide substrate to the cleft and, possibly, locks the enzyme in an “open-jaw”-like conformation making catalysis impossible.

Molecular dynamics study of the recognition of ATP by nucleic acid aptamers

Despite their great success in recognizing small molecules in vitro, nucleic acid aptamers are rarely used in clinical settings. This is partially due to the lack of structure-based mechanistic information. In this work, atomistic molecular dynamics simulations are used to study the static and dynamic supramolecular structures relevant to the process of the wild-type (wt) nucleic acid aptamer recognition and binding of ATP. The effects brought about by mutation of key residues in the recognition site are also explored. The simulations reveal that the aptamer displays a high degree of rigidity and is structurally very little affected by the binding of ATP. Interaction energy decomposition shows that dispersion forces from π-stacking between ATP and the G6 and A23 nucleobases in the aptamer binding site plays a more important role in stabilizing the supramolecular complex, compared to hydrogen-bond interaction between ATP and G22. Moreover, metadynamics simulations show that during the association process, water molecules act as essential bridges connecting ATP with G22, which favors the dynamic stability of the complex. The calculations carried out on three mutated aptamer structures confirm the crucial role of the hydrogen bonds and π-stacking interactions for the binding affinity of the ATP nucleic acid aptamer.

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.

Elucidating the influence of linker histone variants on chromatosome dynamics and energetics

Linker histones are epigenetic regulators that bind to nucleosomes and alter chromatin structures and dynamics. Biophysical studies have revealed two binding modes in the linker histone/nucleosome complex, the chromatosome, where the linker histone is either centered on or askew from the dyad axis. Each has been posited to have distinct effects on chromatin, however the molecular and thermodynamic mechanisms that drive them and their dependence on linker histone compositions remain poorly understood. We present molecular dynamics simulations of chromatosomes with the globular domain of two linker histone variants, generic H1 (genGH1) and H1.0 (GH1.0), to determine how their differences influence chromatosome structures, energetics and dynamics. Results show that both unbound linker histones adopt a single compact conformation. Upon binding, DNA flexibility is reduced, resulting in increased chromatosome compaction. While both variants enthalpically favor on-dyad binding, energetic benefits are significantly higher for GH1.0, suggesting that GH1.0 is more capable than genGH1 of overcoming the large entropic reduction required for on-dyad binding which helps rationalize experiments that have consistently demonstrated GH1.0 in on-dyad states but that show genGH1 in both locations. These simulations highlight the thermodynamic basis for different linker histone binding motifs, and details their physical and chemical effects on chromatosomes.

The shuttling mechanism of foldaxanes: more than just translocation and rotation

Tailoring the structures of nanomachines to achieve specific functions is one of the major challenges in chemistry. Disentangling the different movements of nanomachines is critical to characterize their functions. Here, the motions within one kind of molecular machine, a foldaxane, composed of a foldamer with a spring-like conformation on an axle have been examined at the molecular level. With the aid of molecular dynamics simulations and enhanced sampling methods, the free-energy landscape characterizing the shuttling of the foldaxane has been drawn. The calculated free-energy barrier, amounting to 20.7 kcal mol^-1, is in good agreement with experiments. Further analysis reveals that the predominant contribution to the free-energy barrier stems from the disruption of the hydrogen bonds between the foldamer and the thread. In the absence of hydrogen bonding interactions between the terminals of the foldamer and the thread, shrinkage and swelling movements of the foldamer have been identified and investigated in detail. By deciphering the intricate mechanism of how the foldaxane shuttles, our understanding of motions within molecular machines is expected to be improved, which will, in turn, assist the construction of molecular machines with specific functions.

Molecular Transport across the Ionic Liquid-Aqueous Electrolyte Interface in a MoS2 Nanopore

Nanopore sequencing of DNA has been enabled by the use of a biological enzyme to thread DNA through an engineered biological nanopore while recording the ionic current flowing through the nanopore. Efforts to realize a similar concept using a solid-state nanopore have been met with several technical challenges, one of which is the high speed of DNA translocation and the other the low ionic current contrast among individual nucleotides. A promising avenue to addressing both problems is using an ionic liquid to slow DNA translocation and a tiny nanopore in the MoS₂ membrane to distinguish individual nucleotides. The physical mechanisms enabling these technical advances have remained elusive. Here, we characterize the ion and DNA transport through the ionic liquid/aqueous electrolyte interface, with and without a MoS₂ nanopore, using the all-atom molecular dynamics method. We find that the partial miscibility of the ionic liquid and the aqueous electrolyte considerably alters the physics of the nanopore translocation process. Thus, the interface of the two phases generates a contact potential of 600 mV, the ionic current is dominated by the motion of ionic liquid molecules through the aqueous solution phase, and the DNA nucleotides exhibit preferential partitioning into the aqueous electrolyte, which leads to spontaneous transport of DNA polymers from the ionic liquid to the aqueous solution compartment in the absence of external voltage bias. The complex physics of the two-phase nanopore system offers a multitude of opportunities for extending the functionality of nanopore-sensing platforms.

Cellular absorption of small molecules: free energy landscapes of melatonin binding at phospholipid membranes

Free energy calculations are essential to unveil mechanisms at the atomic scale such as binding of small solutes and their translocation across cell membranes, eventually producing cellular absorption. Melatonin regulates biological rhythms and is directly related to carcinogenesis and neurodegenerative disorders. Free energy landscapes obtained from well-tempered metadynamics simulations precisely describe the characteristics of melatonin binding to specific sites in the membrane and reveal the role of cholesterol in free energy barrier crossing. A specific molecular torsional angle and the distance between melatonin and the center of the membrane along the normal to the membrane Z-axis have been considered as suitable reaction coordinates. Free energy barriers between two particular orientations of the molecular structure (folded and extended) have been found to be of about 18 kJ/mol for z-distances of about 1–2 nm. The ability of cholesterol to expel melatonin out of the internal regions of the membrane towards the interface and the external solvent is explained from a free energy perspective. The calculations reported here offer detailed free energy landscapes of melatonin embedded in model cell membranes and reveal microscopic information on its transition between free energy minima, including the location of relevant transition states, and provide clues on the role of cholesterol in the cellular absorption of small molecules.

Analytical Method Using a Scaled Hypersphere Search for High-Dimensional Metadynamics Simulations

Metadynamics (MTD) is one of the most effective methods for calculating the free energy surface and finding rare events. Nevertheless, numerous studies using MTD have been carried out using 3D or lower dimensional collective variables (CVs), as higher dimensional CVs require costly computational resources and the obtained results are too complex to understand the important events. The latter issue can be conveniently solved by utilizing the free energy reaction network (FERN), which is a graph structure consisting of edges of minimum free energy paths (MFEPs), nodes of equation (EQ) points, and transition state (TS) points. In the present article, a new method for exploring FERNs on high-dimensional CVs using MTD and the scaled hypersphere search (SHS) method is described. A test calculation based on the MTD-SHS simulation of met-enkephalin in explicit water with 7 CVs was conducted. As a result, 889 EQ points and 1805 TS points were found. The MTD-SHS approach can find MFEPs exhaustively; therefore, the FERNs can be estimated without any a priori knowledge of the EQ and TS points.

On the Interpretation of subtilisin Carlsberg Time-Resolved Fluorescence Anisotropy Decays: Modeling with Classical Simulations

In this work, we discuss the challenging time-resolved fluorescence anisotropy of subtilisin Carlsberg (SC), which contains a single Trp residue and is a model fluorescence system. Experimental decay rates and quenching data suggest that the fluorophore should be exposed to water, but the Trp is partially buried in a hydrophobic pocket in the crystallographic structure. In order to study this inconsistency, molecular dynamics simulations were performed to predict the anisotropy decay rates and emission wavelengths of the Trp. We confirmed the inconsistency of the crystallographic structure with the experimentally observed fluorescence data and performed free energy calculations to show that the buried Trp conformation is 2 orders of magnitude (∼3 kcal/mol) more stable than the solvent-exposed one. However, molecular dynamics simulations in which the Trp side chain was restricted to solvent-exposed conformations displayed a maximum Trp emission wavelength shifted toward lower energies and decay rates compatible with the experimentally probed rates. Therefore, if the solvent-exposed conformations are the most important emitters, the experimental anisotropy can be compatibilized with the crystallographic structure. The most likely explanation is that the fluorescence of the most probable conformation in solution, observed in the crystal, is quenched, and this is consistent with the low quantum yield of Trp113 of SC. Additionally, some experiments might have probed denatured or lysed SC structures. SC anisotropy provides an interesting target for the study of fluorescence anisotropy using simulations, which can be used to test and exemplify how modeling can aid the interpretation of experimental data in a system where structure and solution experiments appear to be inconsistent.

Adverse effects of Alport syndrome-related Gly missense mutations on collagen type IV: Insights from molecular simulations and experiments

Patients with Alport syndrome (AS) exhibit blood and elevated protein levels in their urine, inflamed kidneys, and many other abnormalities. AS is attributed to mutations in type IV collagen genes, particularly glycine missense mutations in the collagenous domain of COL4A5 that disrupt common structural motifs in collagen from the repeat (Gly-Xaa-Yaa)_n amino acid sequence. To characterize and elucidate the molecular mechanisms underlying how AS-related mutations perturb the structure and function of type IV collagen, experimental studies and molecular simulations were integrated to investigate the structure, stability, protease sensitivity, and integrin binding affinity of collagen-like proteins containing amino acid sequences from the α5(IV) chain and AS-related Gly missense mutations. We show adverse effects where (i) three AS-related Gly missense mutations significantly reduced the structural stability of the collagen in terms of decreased melting temperatures and calorimetric enthalpies, in conjunction with a collective drop in the external work needed to unfold the peptides containing mutation sequences; (ii) due to local unwinding around the sites of mutations, these triple helical peptides were also degraded more rapidly by trypsin and chymotrypsin, as these enzymes could access the collagenous triple helix more easily and increase the number of contacts; (iii) the mutations further abolished the ability of the recombinant collagens to bind to integrins and greatly reduced the binding affinities between collagen and integrins, thus preventing cells from adhering to these mutants. Our unified experimental and computational approach provided underlying insights needed to guide potential therapies for AS that ameliorate the adverse effects from AS disease onset and progression.

Capturing Structural Transitions in Surfactant Adsorption Isotherms at Solid/Solution Interfaces

Although adsorption isotherms of surfactants are critical in determining the relationship between interfacial properties and structures of surfactants, providing quantitative predictions of the isotherms remains challenging. This is especially true for adsorption at hard interfaces such as on two-dimensional (2D) layered materials or on nanoparticles where simulation techniques developed for fluid-fluid interfaces that dynamically change surface properties by adjusting unit cells do not apply. In this work, we predict nonideal adsorption at a solid-solution interface with a molecular thermodynamic theory (MTT) model that utilizes molecular dynamics (MD) simulations for the determination of free-energy parameters in the MTT. Furthermore, the MD/MTT model provides atomistic insights into the nonideal behavior of surfactants by capturing structural phases of the surfactants at the interface. Our approach captures structural transitions from the ideal state at low concentrations and then to the critical surface aggregation concentration (CSAC) and finally through the critical micelle concentration (CMC). We validate our model against the original MTT model by comparing predicted adsorption isotherms of a simplified surfactant system from both approaches. We further substantiate the applicability of our model in complex systems by providing adsorption isotherms in an aqueous sodium dodecyl sulfate (SDS)-graphene system, in good agreement with the experimental observations of the CSAC for the same system.

Atomistic Simulations of COSAN: Amphiphiles without a Head-and-Tail Design Display "Head and Tail" Surfactant Behavior

Cobaltabisdicarbollide (COSAN) anions have an unexpectedly rich self-assembly behavior, which can lead to vesicles and micelles without having a classical surfactant molecular architecture. This was rationalized by the introduction of new terminology and novel driving forces. A key aspect in the interpretation of COSAN behavior is the assumption that the most stable form of these ions is the transoid rotamer, which lacks a "hydrophilic head" and a "hydrophobic tail". Using implicit solvent DFT calculations and MD simulations we show that in water, 1) the cisoid rotamer is the most stable form of COSAN and 2) this cisoid rotamer has a well-defined hydrophilic polar region ("head") and a hydrophobic apolar region ("tail"). In addition, our simulations show that the properties of this rotamer in water (interfacial affinity, micellization) match those expected for a classical surfactant. Therefore, we conclude that the experimental results for the COSAN ions can now be understood in terms of its amphiphilic molecular architecture.

Cucurbit[8]uril-mediated pseudo[2,3]rotaxanes

Pseudo[2,3]rotaxanes have been successfully fabricated by the complexation of cucurbit[8]uril (CB[8]) macrocycles with extended viologen derivatives. Two design rules enable the incorporation of a third CB[8] onto a recently reported pseudo[2,2]rotaxane. Incorporation of a third macrocycle confines the dimeric stacking of chromophores into specific alignment, leading to effective electron-delocalisation along their long molecular axis.

Synergistic effect of supercritical CO2 and organic solvent on exfoliation of graphene: experiment and atomistic simulation studies

In this work, experiments and molecular dynamics (MD) simulations are carried out to explore the synergistic effect of supercritical CO2 (scCO2) and organic solvent on intercalation and exfoliation of graphene. Experimental characterizations via transmission electron microscopy, atomic force microscopy and Raman spectroscopy indicate that by combining scCO2 and organic solvent (N-methylpyrrolidone, NMP), few-layer graphene is successfully exfoliated from graphite, among which over 30% is 1-4 layers, and 55% is 5-8 layers. Systematic experiments have shown that compared with pure scCO2 or NMP, the mixed scCO2 and NMP can significantly increase the amount of graphene and the rate of few-layer graphene, and the optimum volume fraction of NMP is 25%. Parallel MD simulations indicate that the scCO2 molecules first diffuse into the interlayer of graphite, and then the larger NMP molecules insert as wedges and further expand interlayer spacing, promoting intercalation and exfoliation. The iteration of scCO2 diffusion and the NMP wedge can generate positive feedback to improve the exfoliation productivity and efficiency. This work explores the synergistic effect of scCO2 and NMP on the exfoliation of graphene, which may provide useful insights for exfoliation of other two dimensional materials.

Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning

Computing equilibrium states in condensed-matter many-body systems, such as solvated proteins, is a long-standing challenge. Lacking methods for generating statistically independent equilibrium samples in “one shot,” vast computational effort is invested for simulating these systems in small steps, e.g., using molecular dynamics. Combining deep learning and statistical mechanics, we developed Boltzmann generators, which are shown to generate unbiased one-shot equilibrium samples of representative condensed-matter systems and proteins. Boltzmann generators use neural networks to learn a coordinate transformation of the complex configurational equilibrium distribution to a distribution that can be easily sampled. Accurate computation of free-energy differences and discovery of new configurations are demonstrated, providing a statistical mechanics tool that can avoid rare events during sampling without prior knowledge of reaction coordinates.

Modulation of Antimicrobial Peptide Potency in Stressed Lipid Bilayers

It is shown that the tendency of an archetypal antimicrobial peptide to insert into and perforate a simple lipid bilayer is strongly modulated by tensile stress in the membrane. The results, obtained through molecular dynamics simulations, have been demonstrated with several lipid compositions and appear to be general, although quantitative details differ. The findings imply that the potency of antimicrobial peptides may not be a purely intrinsic chemical property and, instead, depends on the mechanical state of the target membrane.

Molecular insight into the wetting behavior and amphiphilic character of cellulose nanocrystals

The study of nanocellulose is a field of growing interest due to its many applications and its use in the development of biocompatible and eco-friendly materials. In spite of the vast number of studies in the field, many questions about the role of the molecular structure in the properties of cellulose are still subject of debate. One of these fundamental questions is the possible amphiphilic nature of cellulose and the relative role of hydrogen bonding and hydrophobic effect on the interactions of cellulose. In this work we present an extensive molecular dynamics simulation study of this question by analyzing the wetting of cellulose with water and organic solvent, its interaction with hydrophilic and hydrophobic ions and its interaction with a protein (human epidermal growth factor, hEGF). We consider two characteristic cellulose crystal planes of Iβ cellulose with very different roughness, different hydrogen bonding capability and different exposure of cellulose hydrophobic groups (the (010) plane which has exposed -OH groups and the (100) plane with buried -OH groups). Our results show that both surfaces are simultaneously hydrophilic and lipophilic, with both surfaces having very similar contact angles. In spite of the global similarity of wetting of both surfaces, the molecular details of wetting are very different and substantial local wetting heterogeneities (which strongly depend on the surface) appear for both solvents. We also observe a weak interaction of both surfaces with hydrophobic and hydrophilic solutes. These weak interactions are attributed to the simultaneous lipophilic and hydrophilic character of both (100) and (010) cellulose surfaces. Interestingly, we found a substantial interaction of both cellulose planes with polar and apolar residues of the hEGF protein.

Q817G mutation in phosphodiesterase type 5: Conformational analysis and dissociation profile of the inhibitor Tadalafil

Phosphodiesterase type 5 (PDE-5) is an important enzyme involved in the hydrolysis of cyclic guanosine monophosphate (cGMP) to guanosine monophosphate (GMP). The inhibition of this protein leads to the accumulation of cGMP in cells with various biological and therapeutic effects. Several PDE-5 inhibitors exist, with Tadalafil being one of the most commonly studied and used in clinical therapy. In this study, we applied Molecular Dynamics simulations coupled to the ABF (Adaptive Biasing Force) method to study the effect of the mutation on the Gln817 residue (Q817G). The results of the free energy profiles made clear that the affinity of the inhibitor for PDE-5 is dependent on the amino acid residue Gln817. The hydrogen bond made between the side chain of glutamine and the indole ring of Tadalafil results in the stabilization of the ligand in the catalytic site. Despite the prominent role of this interaction, it is important to highlight the contribution of other residues of the catalytic domain for the stabilization of the compound, due to the set of polar, hydrophobic and electrostatic interactions performed by specific amino acid residues.