Optimized reaction coordinates for analysis of enhanced sampling.

Atomistic simulations of biological processes offer insights at a high level of spatial and temporal resolution, but accelerated sampling is often required for probing timescales of biologically relevant processes. The resulting data need to be statistically reweighted and condensed in a concise yet faithful manner to facilitate interpretation. Here, we provide evidence that a recently proposed approach for the unsupervised determination of optimized reaction coordinate (RC) can be used for both analysis and reweighting of such data. We first show that for a peptide interconverting between helical and collapsed configurations, the optimal RC permits efficient reconstruction of equilibrium properties from enhanced sampling trajectories. Upon RC-reweighting, kinetic rate constants and free energy profiles are in good agreement with values obtained from equilibrium simulations. In a more challenging test, we apply the method to enhanced sampling simulations of the unbinding of an acetylated lysine-containing tripeptide from the bromodomain of ATAD2. The complexity of this system allows us to investigate the strengths and limitations of these RCs. Overall, the findings presented here underline the potential of the unsupervised determination of reaction coordinates and the synergy with orthogonal analysis methods, such as Markov state models and SAPPHIRE analysis.

Antibody binding modulates the dynamics of the membrane-bound prion protein.

Misfolding of the cellular prion protein (PrPC) is associated with lethal neurodegeneration. PrPC consists of a flexible tail (residues 23-123) and a globular domain (residues 124-231) whose C-terminal end is anchored to the cell membrane. The neurotoxic antibody POM1 and the innocuous antibody POM6 recognize the globular domain. Experimental evidence indicates that POM1 binding to PrPC emulates the influence on PrPC of the misfolded prion protein (PrPSc) while the binding of POM6 has the opposite biological response. Little is known about the potential interactions between flexible tail, globular domain, and the membrane. Here, we used atomistic simulations to investigate how these interactions are modulated by the binding of the Fab fragments of POM1 and POM6 to PrPC and by interstitial sequence truncations to the flexible tail. The simulations show that the binding of the antibodies restricts the range of orientations of the globular domain with respect to the membrane and decreases the distance between tail and membrane. Five of the six sequence truncations influence only marginally this distance and the contact patterns between tail and globular domain. The only exception is a truncation coupled to a charge inversion mutation of four N-terminal residues, which increases the distance of the flexible tail from the membrane. The interactions of the flexible tail and globular domain are modulated differently by the two antibodies.

An ABSINTH-Based Protocol for Predicting Binding Affinities between Proteins and Small Molecules.

The core task in computational drug discovery is to accurately predict binding free energies in receptor-ligand systems for large libraries of putative binders. Here, the ABSINTH implicit solvent model and force field are extended to describe small, organic molecules and their interactions with proteins. We show that an automatic pipeline based on partitioning arbitrary molecules into substructures corresponding to model compounds with known free energies of solvation can be combined with the CHARMM general force field into a method that is successful at the two important challenges a scoring function faces in virtual screening work flows: it ranks known binders with correlation values rivaling that of comparable state-of-the-art methods and it enriches true binders in a set of decoys. Our protocol introduces innovative modifications to common virtual screening workflows, notably the use of explicit ions as competitors and the integration over multiple protein and ligand species differing in their protonation states. We demonstrate the value of modifications to both the protocol and ABSINTH itself. We conclude by discussing the limitations of high-throughput implicit methods such as the one proposed here.

On the removal of initial state bias from simulation data.

Classical atomistic simulations of biomolecules play an increasingly important role in molecular life science. The structure of current computing architectures favors methods that run multiple trajectories at once without requiring extensive communication between them. Many advanced sampling strategies in the field fit this mold. These approaches often rely on an adaptive logic and create ensembles of comparatively short trajectories whose starting points are not distributed according to the correct Boltzmann weights. This type of bias is notoriously difficult to remove, and Markov state models (MSMs) are one of the few strategies available for recovering the correct kinetics and thermodynamics from these ensembles of trajectories. In this contribution, we analyze the performance of MSMs in the thermodynamic reweighting task for a hierarchical set of systems. We show that MSMs can be rigorous tools to recover the correct equilibrium distribution for systems of sufficiently low dimensionality. This is conditional upon not tampering with local flux imbalances found in the data. For a real-world application, we find that a pure likelihood-based inference of the transition matrix produces the best results. The removal of the bias is incomplete, however, and for this system, all tested MSMs are outperformed by an alternative albeit less general approach rooted in the ideas of statistical resampling. We conclude by formulating some recommendations for how to address the reweighting issue in practice.

The ATAD2 bromodomain binds different acetylation marks on the histone H4 in similar fuzzy complexes.

Bromodomains are protein modules adopting conserved helix bundle folds. Some bromodomain-containing proteins, such as ATPase family AAA domain-containing protein 2 (ATAD2), isoform A, have attracted much interest because they are overexpressed in many types of cancer. Bromodomains bind to acetylated lysine residues on histone tails and thereby facilitate the reading of the histone code. Epigenetic regulators in general have been implicated as indicators, mediators, or causes of a large number of diseases and disorders. To interfere with or modulate these processes, it is therefore of fundamental interest to understand the molecular mechanisms by which epigenetic regulation occurs. Here, we present results from molecular dynamics simulations of a doubly acetylated histone H4 peptide bound to the bromodomain of ATAD2 (hereafter referred to as ATAD2A). These simulations revealed how the flexibility of ATAD2A's major loop, the so-called ZA loop, creates an adaptable interface that preserves the disorder of both peptide and loop in the bound state. We further demonstrate that the binding involves an almost identical average pattern of interactions irrespective of which acetyl mark is inserted into the pocket. In conjunction with a likely mechanism of electrostatically driven recruitment, our simulation results highlight how the bromodomain is built toward promiscuous binding with low specificity. In conclusion, the simulations indicate that disorder and electrostatic steering function jointly to recruit ATAD2A to the histone core and that these fuzzy interactions may promote cooperativity between nearby epigenetic marks.

Precise estimation of transfer free energies for ionic species between similar media.

Transfer properties and partition coefficients for individual ions are relevant in a variety of scientific and engineering contexts, such as predicting the effects of different electrolytes on biomacromolecules in a preferential interaction sense or predicting the distribution of heavy metal ions in soils, rivers, etc. Computer simulations allow free energies of transfer to be estimated by considering single ions explicitly. When the two media under consideration are similar to each other regarding ion solvation, the resultant free energies are small in absolute magnitude. In these cases, it is advisable to simulate the transfer process directly. Here, we demonstrate how this can be achieved using two-dimensional umbrella sampling in conjunction with canonical ensemble molecular dynamics simulations where two liquid media are in direct contact. By calculating full two-dimensional potentials of mean force, these simulations allow the estimation of single-ion transfer free energies by integrating this surface accordingly. We report statistical accuracies to highlight that very high precision is achieved and needed to make even just qualitative statements about the transfer process. We close by discussing implications of our results for the specific case considered: the transfer of polypeptide side chain analogs from water to aqueous denaturant solutions.

Focused conformational sampling in proteins.

A detailed understanding of the conformational dynamics of biological molecules is difficult to obtain by experimental techniques due to resolution limitations in both time and space. Computer simulations avoid these in theory but are often too short to sample rare events reliably. Here we show that the progress index-guided sampling (PIGS) protocol can be used to enhance the sampling of rare events in selected parts of biomolecules without perturbing the remainder of the system. The method is very easy to use as it only requires as essential input a set of several features representing the parts of interest sufficiently. In this feature space, new states are discovered by spontaneous fluctuations alone and in unsupervised fashion. Because there are no energetic biases acting on phase space variables or projections thereof, the trajectories PIGS generates can be analyzed directly in the framework of transition networks. We demonstrate the possibility and usefulness of such focused explorations of biomolecules with two loops that are part of the binding sites of bromodomains, a family of epigenetic "reader" modules. This real-life application uncovers states that are structurally and kinetically far away from the initial crystallographic structures and are also metastable. Representative conformations are intended to be used in future high-throughput virtual screening campaigns.

A molecular simulation protocol to avoid sampling redundancy and discover new states.

BACKGROUND: For biomacromolecules or their assemblies, experimental knowledge is often restricted to specific states. Ambiguity pervades simulations of these complex systems because there is no prior knowledge of relevant phase space domains, and sampling recurrence is difficult to achieve. In molecular dynamics methods, ruggedness of the free energy surface exacerbates this problem by slowing down the unbiased exploration of phase space. Sampling is inefficient if dwell times in metastable states are large. METHODS: We suggest a heuristic algorithm to terminate and reseed trajectories run in multiple copies in parallel. It uses a recent method to order snapshots, which provides notions of "interesting" and "unique" for individual simulations. We define criteria to guide the reseeding of runs from more "interesting" points if they sample overlapping regions of phase space. RESULTS: Using a pedagogical example and an α-helical peptide, the approach is demonstrated to amplify the rate of exploration of phase space and to discover metastable states not found by conventional sampling schemes. Evidence is provided that accurate kinetics and pathways can be extracted from the simulations. CONCLUSIONS: The method, termed PIGS for Progress Index Guided Sampling, proceeds in unsupervised fashion, is scalable, and benefits synergistically from larger numbers of replicas. Results confirm that the underlying ideas are appropriate and sufficient to enhance sampling. GENERAL SIGNIFICANCE: In molecular simulations, errors caused by not exploring relevant domains in phase space are always unquantifiable and can be arbitrarily large. Our protocol adds to the toolkit available to researchers in reducing these types of errors. This article is part of a Special Issue entitled "Recent developments of molecular dynamics".

Quantitative assessments of the distinct contributions of polypeptide backbone amides versus side chain groups to chain expansion via chemical denaturation.

In aqueous solutions with high concentrations of chemical denaturants such as urea and guanidinium chloride (GdmCl) proteins expand to populate heterogeneous conformational ensembles. These denaturing environments are thought to be good solvents for generic protein sequences because properties of conformational distributions align with those of canonical random coils. Previous studies showed that water is a poor solvent for polypeptide backbones, and therefore, backbones form collapsed globular structures in aqueous solvents. Here, we ask if polypeptide backbones can intrinsically undergo the requisite chain expansion in aqueous solutions with high concentrations of urea and GdmCl. We answer this question using a combination of molecular dynamics simulations and fluorescence correlation spectroscopy. We find that the degree of backbone expansion is minimal in aqueous solutions with high concentrations of denaturants. Instead, polypeptide backbones sample conformations that are denaturant-specific mixtures of coils and globules, with a persistent preference for globules. Therefore, typical denaturing environments cannot be classified as good solvents for polypeptide backbones. How then do generic protein sequences expand in denaturing environments? To answer this question, we investigated the effects of side chains using simulations of two archetypal sequences with amino acid compositions that are mixtures of charged, hydrophobic, and polar groups. We find that side chains lower the effective concentration of backbone amides in water leading to an intrinsic expansion of polypeptide backbones in the absence of denaturants. Additional dilution of the effective concentration of backbone amides is achieved through preferential interactions with denaturants. These effects lead to conformational statistics in denaturing environments that are congruent with those of canonical random coils. Our results highlight the role of side chain-mediated interactions as determinants of the conformational properties of unfolded states in water and in influencing chain expansion upon denaturation.

A simple molecular mechanics integrator in mixed rigid body and dihedral angle space.

We propose a numerical scheme to integrate equations of motion in a mixed space of rigid-body and dihedral angle coordinates. The focus of the presentation is biomolecular systems and the framework is applicable to polymers with tree-like topology. By approximating the effective mass matrix as diagonal and lumping all bias torques into the time dependencies of the diagonal elements, we take advantage of the formal decoupling of individual equations of motion. We impose energy conservation independently for every degree of freedom and this is used to derive a numerical integration scheme. The cost of all auxiliary operations is linear in the number of atoms. By coupling the scheme to one of two popular thermostats, we extend the method to sample constant temperature ensembles. We demonstrate that the integrator of choice yields satisfactory stability and is free of mass-metric tensor artifacts, which is expected by construction of the algorithm. Two fundamentally different systems, viz., liquid water and an α-helical peptide in a continuum solvent are used to establish the applicability of our method to a wide range of problems. The resultant constant temperature ensembles are shown to be thermodynamically accurate. The latter relies on detailed, quantitative comparisons to data from reference sampling schemes operating on exactly the same sets of degrees of freedom.

Net charge per residue modulates conformational ensembles of intrinsically disordered proteins.

Intrinsically disordered proteins (IDPs) adopt heterogeneous ensembles of conformations under physiological conditions. Understanding the relationship between amino acid sequence and conformational ensembles of IDPs can help clarify the role of disorder in physiological function. Recent studies revealed that polar IDPs favor collapsed ensembles in water despite the absence of hydrophobic groups--a result that holds for polypeptide backbones as well. By studying highly charged polypeptides, a different archetype of IDPs, we assess how charge content modulates the intrinsic preference of polypeptide backbones for collapsed structures. We characterized conformational ensembles for a set of protamines in aqueous milieus using molecular simulations and fluorescence measurements. Protamines are arginine-rich IDPs involved in the condensation of chromatin during spermatogenesis. Simulations based on the ABSINTH implicit solvation model predict the existence of a globule-to-coil transition, with net charge per residue serving as the discriminating order parameter. The transition is supported by quantitative agreement between simulation and experiment. Local conformational preferences partially explain the observed trends of polymeric properties. Our results lead to the proposal of a schematic protein phase diagram that should enable prediction of polymeric attributes for IDP conformational ensembles using easily calculated physicochemical properties of amino acid sequences. Although sequence composition allows the prediction of polymeric properties, interresidue contact preferences of protamines with similar polymeric attributes suggest that certain details of conformational ensembles depend on the sequence. This provides a plausible mechanism for specificity in the functions of IDPs.

Disordered binding of small molecules to Aß(12-28).

In recent years, an increasing number of small molecules and short peptides have been identified that interfere with aggregation and/or oligomerization of the Alzheimer ß-amyloid peptide (Aß). Many of them possess aromatic moieties, suggesting a dominant role for those in interacting with Aß along various stages of the aggregation process. In this study, we attempt to elucidate whether interactions of such aromatic inhibitors with monomeric Aß(12-28) point to a common mechanism of action by performing atomistic molecular dynamics simulations at equilibrium. Our results suggest that, independently of the presence of inhibitors, monomeric Aß(12-28) populates a partially collapsed ensemble that is largely devoid of canonical secondary structure at 300 K and neutral pH. The small molecules have different affinities for Aß(12-28) that can be partially rationalized by the balance of aromatic and charged moieties constituting the molecules. There are no predominant binding modes, although aggregation inhibitors preferentially interact with the N-terminal portion of the fragment (residues 13-20). Analysis of the free energy landscape of Aß(12-28) reveals differences highlighted by altered populations of a looplike conformer in the presence of inhibitors. We conclude that intrinsic disorder of Aß persists at the level of binding small molecules and that inhibitors can significantly alter properties of monomeric Aß via multiple routes of differing specificity.

Spiking burstiness and working memory in the human medial temporal lobe.

Persistent activity has commonly been considered to be a hallmark of working memory (WM). Recent evidence indicates that neuronal discharges in the medial temporal lobe (MTL) are compatible with WM neural patterns observed in cortical areas. However, the characterization of this activity rarely consists of measurements other than firing rates of single neurons. Moreover, a varied repertoire of firing dynamics has been reported in the MTL regions, which motivate the more detailed examination of the relationships between WM processes and discharge patterns undertaken here. Specifically, we investigate' at different resolution levels, firing irregularities in electrode recordings from the hippocampus, amygdala, and the entorhinal cortex of epileptic patients during a WM task. We show that some types of (ir)regularities predict response times of the patients depending on the trial periods under consideration. Prominent burst activity at the population level is observed in the amygdala and entorhinal cortex during memory retrieval. In general, regular and bursty neurons contribute to the decoding of the memory load, yet they display important differences across the three anatomical areas. Our results suggest that nonrandom (non-Poisson) patterns are relevant for WM, which calls for the development and use of statistics complementary to mere spike counts.

Unsupervised Methods for Detection of Neural States: Case Study of Hippocampal-Amygdala Interactions.

The hippocampus and amygdala are functionally coupled brain regions that play a crucial role in processes involving memory and learning. Because interareal communication has been reported both during specific sleep stages and in awake, behaving animals, these brain regions can serve as an archetype to establish that measuring functional interactions is important for comprehending neural systems. To this end, we analyze here a public dataset of local field potentials (LFPs) recorded in rats simultaneously from the hippocampus and amygdala during different behaviors. Employing a specific, time-lagged embedding technique, named topological causality (TC), we infer directed interactions between the LFP band powers of the two regions across six frequency bands in a time-resolved manner. The combined power and interaction signals are processed with our own unsupervised tools developed originally for the analysis of molecular dynamics simulations to effectively visualize and identify putative, neural states that are visited by the animals repeatedly. Our proposed methodology minimizes impositions onto the data, such as isolating specific epochs, or averaging across externally annotated behavioral stages, and succeeds in separating internal states by external labels such as sleep or stimulus events. We show that this works better for two of the three rats we analyzed, and highlight the need to acknowledge individuality in analyses of this type. Importantly, we demonstrate that the quantification of functional interactions is a significant factor in discriminating these external labels, and we suggest our methodology as a general tool for large, multisite recordings.

Sapphire-Based Clustering.

Molecular dynamics simulations are a popular means to study biomolecules, but it is often difficult to gain insights from the trajectories due to their large size, in both time and number of features. The Sapphire (States And Pathways Projected with HIgh REsolution) plot allows a direct visual inference of the dominant states visited by high-dimensional systems and how they are interconnected in time. Here, we extend this visual inference into a clustering algorithm. Specifically, the automatic procedure derives from the Sapphire plot states that are kinetically homogeneous, structurally annotated, and of tunable granularity. We provide a relative assessment of the kinetic fidelity of the Sapphire-based partitioning in comparison to popular clustering methods. This assessment is carried out on trajectories of n-butane, a ß-sheet peptide, and the small protein BPTI. We conclude with an application of our approach to a recent 100 µs trajectory of the main protease of SARS-CoV-2.

Thermodynamics of beta-sheet formation in polyglutamine.

The role of beta-sheets in the early stages of protein aggregation, specifically amyloid formation, remains unclear. Interpretations of kinetic data have led to a specific model for the role of beta-sheets in polyglutamine aggregation. According to this model, monomeric polyglutamine, which is intrinsically disordered, goes through a rare conversion into an ordered, metastable, beta-sheeted state that nucleates aggregation. It has also been proposed that the probability of forming the critical nucleus, a specific beta-sheet conformation for the monomer, increases with increasing chain length. Here, we test this model using molecular simulations. We quantified free energy profiles in terms of beta-content for monomeric polyglutamine as a function of chain length. In accord with estimates from experimental data, the free energy penalties for forming beta-rich states are in the 10-20 kcal/mol range. However, the length dependence of these free energy penalties does not mirror interpretations of kinetic data. In addition, although homodimerization of disordered molecules is spontaneous, the imposition of conformational restraints on polyglutamine molecules does not enhance the spontaneity of intermolecular associations. Our data lead to the proposal that beta-sheet formation is an attribute of peptide-rich phases such as high molecular weight aggregates rather than monomers or oligomers.

ABSINTH: a new continuum solvation model for simulations of polypeptides in aqueous solutions.

A new implicit solvation model for use in Monte Carlo simulations of polypeptides is introduced. The model is termed ABSINTH for self-Assembly of Biomolecules Studied by an Implicit, Novel, and Tunable Hamiltonian. It is designed primarily for simulating conformational equilibria and oligomerization reactions of intrinsically disordered proteins in aqueous solutions. The paradigm for ABSINTH is conceptually similar to the EEF1 model of Lazaridis and Karplus (Proteins 1999, 35, 133). In ABSINTH, the transfer of a polypeptide solute from the gas phase into a continuum solvent is the sum of a direct mean field interaction (DMFI), and a term to model the screening of polar interactions. Polypeptide solutes are decomposed into a set of distinct solvation groups. The DMFI is a sum of contributions from each of the solvation groups, which are analogs of model compounds. Continuum-mediated screening of electrostatic interactions is achieved using a framework similar to the one used for the DMFI. Promising results are shown for a set of test cases. These include the calculation of NMR coupling constants for short peptides, the assessment of the thermal stability of two small proteins, reversible folding of both an alpha-helix and a beta-hairpin forming peptide, and the polymeric properties of intrinsically disordered polyglutamine peptides of varying lengths. The tests reveal that the computational expense for simulations with the ABSINTH implicit solvation model increase by a factor that is in the range of 2.5-5.0 with respect to gas-phase calculations.

Unsupervised identification of states from voltage recordings of neural networks.

BACKGROUND: Modern techniques for multi-neuronal recording produce large amounts of data. There is no automatic procedure for the identification of states in recurrent voltage patterns. NEW METHOD: We propose NetSAP (Network States And Pathways), a data-driven analysis method that is able to recognize multi-neuron voltage patterns (states). To capture the subtle differences between snapshots in voltage recordings, NetSAP infers the underlying functional neural network in a time-resolved manner with a sliding window approach. Then NetSAP identifies states from a reordering of the time series of inferred networks according to a user-defined metric. The procedure for unsupervised identification of states was developed originally for the analysis of molecular dynamics simulations of proteins. RESULTS: We tested NetSAP on neural network simulations of GABAergic inhibitory interneurons. Most simulation parameters are chosen to reproduce literature observations, and we keep noise terms as control parameters to regulate the coherence of the simulated signals. NetSAP is able to identify multiple states even in the case of high internal noise and low signal coherence. We provide evidence that NetSAP is robust for networks with up to about 50% of the neurons spiking randomly. NetSAP is scalable and its code is open source. COMPARISON WITH EXISTING METHODS: NetSAP outperforms common analysis techniques, such as PCA and k-means clustering, on a simulated recording of voltage traces of 50 neurons. CONCLUSIONS: NetSAP analysis is an efficient tool to identify voltage patterns from neuronal recordings.

A polymer physics perspective on driving forces and mechanisms for protein aggregation.

Protein aggregation is a commonly occurring problem in biology. Cells have evolved stress-response mechanisms to cope with problems posed by protein aggregation. Yet, these quality control mechanisms are overwhelmed by chronic aggregation-related stress and the resultant consequences of aggregation become toxic to cells. As a result, a variety of systemic and neurodegenerative diseases are associated with various aspects of protein aggregation and rational approaches to either inhibit aggregation or manipulate the pathways to aggregation might lead to an alleviation of disease phenotypes. To develop such approaches, one needs a rigorous and quantitative understanding of protein aggregation. Much work has been done in this area. However, several unanswered questions linger, and these pertain primarily to the actual mechanism of aggregation as well as to the types of inter-molecular associations and intramolecular fluctuations realized at low protein concentrations. It has been suggested that the concepts underlying protein aggregation are similar to those used to describe the aggregation of synthetic polymers. Following this suggestion, the relevant concepts of polymer aggregation are introduced. The focus is on explaining the driving forces for polymer aggregation and how these driving forces vary with chain length and solution conditions. It is widely accepted that protein aggregation is a nucleation-dependent process. This view is based mainly on the presence of long times for the accumulation of aggregates and the elimination of these lag times with "seeds". In this sense, protein aggregation is viewed as being analogous to the aggregation of colloidal particles. The theories for polymer aggregation reviewed in this work suggest an alternative mechanism for the origin of long lag times in protein aggregation. The proposed mechanism derives from the recognition that polymers have unique dynamics that distinguish them from other aggregation-prone systems such as colloidal particles.