Building quantum mechanics quality force fields of proteins with the generalized energy-based fragmentation approach and machine learning.

We combined our generalized energy-based fragmentation (GEBF) approach and machine learning (ML) technique to construct quantum mechanics (QM) quality force fields for proteins. In our scheme, the training sets for a protein are only constructed from its small subsystems, which capture all short-range interactions in the target system. The energy of a given protein is expressed as the summation of atomic contributions from QM calculations of various subsystems, corrected by long-range Coulomb and van der Waals interactions. With the Gaussian approximation potential (GAP) method, our protocol can automatically generate training sets with high efficiency. To facilitate the construction of training sets for proteins, we store all trained subsystem data in a library. If subsystems in the library are detected in a new protein, corresponding datasets can be directly reused as a part of the training set on this new protein. With two polypeptides, 4ZNN and 1XQ8 segment, as examples, the energies and forces predicted by GEBF-GAP are in good agreement with those from conventional QM calculations, and dihedral angle distributions from GEBF-GAP molecular dynamics (MD) simulations can also well reproduce those from ab initio MD simulations. In addition, with the training set generated from GEBF-GAP, we also demonstrate that GEBF-ML force fields constructed by neural network (NN) methods can also show QM quality. Therefore, the present work provides an efficient and systematic way to build QM quality force fields for biological systems.

Markov state models from hierarchical density-based assignment.

Markov state models (MSMs) have become one of the preferred methods for the analysis and interpretation of molecular dynamics (MD) simulations of conformational transitions in biopolymers. While there is great variation in terms of implementation, a well-defined workflow involving multiple steps is often adopted. Typically, molecular coordinates are first subjected to dimensionality reduction and then clustered into small "microstates," which are subsequently lumped into "macrostates" using the information from the slowest eigenmodes. However, the microstate dynamics is often non-Markovian, and long lag times are required to converge the relevant slow dynamics in the MSM. Here, we propose a variation on this typical workflow, taking advantage of hierarchical density-based clustering. When applied to simulation data, this type of clustering separates high population regions of conformational space from others that are rarely visited. In this way, density-based clustering naturally implements assignment of the data based on transitions between metastable states, resulting in a core-set MSM. As a result, the state definition becomes more consistent with the assumption of Markovianity, and the timescales of the slow dynamics of the system are recovered more effectively. We present results of this simplified workflow for a model potential and MD simulations of the alanine dipeptide and the FiP35 WW domain.

The synergistic mechanisms of apo-ferritin structural transitions and Au(iii) ion transportation: molecular dynamics simulations with the Markov state model.

Due to its unique structure, recent years have witnessed the use of apo-ferritin to accumulate various non-natural metal ions as a scaffold for nanomaterial synthesis. However, the transport mechanism of metal ions into the cavity of apo-ferritin is still unclear, limiting the rational design and controllable preparation of nanomaterials. Here, we conducted all-atom classical molecular dynamics (MD) simulations combined with Markov state models (MSMs) to explore the transportation behavior of Au(iii) ions. We exhibited the complete transportation paths of Au(iii) from solution into the apo-ferritin cage at the atomic level. We also revealed that the transportation of Au(iii) ions is accompanied by coupled protein structural changes. It is shown that the 3-fold axis channel serves as the only entrance with the longest residence time of Au(iii) ions. Besides, there are eight binding clusters and five 3-fold structural metastable states, which are important during Au(iii) transportation. The conformational changes of His118, Asp127, and Glu130, acting as doors, were observed to highly correlate with the Au(iii) ion's position. The MSM analysis and Potential Mean Force (PMF) calculation suggest a remarkable energy barrier near Glu130, making it the rate-limiting step of the whole process. The dominant transportation pathway is from cluster 3 in the 3-fold channel to the inner cavity to cluster 5 on the inner surface, and then to cluster 6. These findings provide inspiration and theoretical guidance for the further rational design and preparation of new nanomaterials using apo-ferritin.

Methodology for rigorous modeling of protein conformational changes by Rosetta using DEER distance restraints.

We describe an approach for integrating distance restraints from Double Electron-Electron Resonance (DEER) spectroscopy into Rosetta with the purpose of modeling alternative protein conformations from an initial experimental structure. Fundamental to this approach is a multilateration algorithm that harnesses sets of interconnected spin label pairs to identify optimal rotamer ensembles at each residue that fit the DEER decay in the time domain. Benchmarked relative to data analysis packages, the algorithm yields comparable distance distributions with the advantage that fitting the DEER decay and rotamer ensemble optimization are coupled. We demonstrate this approach by modeling the protonation-dependent transition of the multidrug transporter PfMATE to an inward facing conformation with a deviation to the experimental structure of less than 2Å Cα RMSD. By decreasing spin label rotamer entropy, this approach engenders more accurate Rosetta models that are also more closely clustered, thus setting the stage for more robust modeling of protein conformational changes.

Exploring ligand binding pathways on proteins using hypersound-accelerated molecular dynamics.

Capturing the dynamic processes of biomolecular systems in atomistic detail remains difficult despite recent experimental advances. Although molecular dynamics (MD) techniques enable atomic-level observations, simulations of "slow" biomolecular processes (with timescales longer than submilliseconds) are challenging because of current computer speed limitations. Therefore, we developed a method to accelerate MD simulations by high-frequency ultrasound perturbation. The binding events between the protein CDK2 and its small-molecule inhibitors were nearly undetectable in 100-ns conventional MD, but the method successfully accelerated their slow binding rates by up to 10-20 times. Hypersound-accelerated MD simulations revealed a variety of microscopic kinetic features of the inhibitors on the protein surface, such as the existence of different binding pathways to the active site. Moreover, the simulations allowed the estimation of the corresponding kinetic parameters and exploring other druggable pockets. This method can thus provide deeper insight into the microscopic interactions controlling biomolecular processes.

Comprehensive characterization of oligosaccharide conformational ensembles with conformer classification by free-energy landscape via reproductive kernel Hilbert space.

Oligosaccharides play versatile roles in various biological systems but are difficult to characterize from a structural viewpoint due to their remarkable degrees of freedom in internal motion. Therefore, molecular dynamics simulations have been widely used to delineate the dynamic conformations of oligosaccharides. However, hardly any methods have thus far been available for the comprehensive characterization of simulation-derived conformational ensembles of oligosaccharides. In this research, we attempted to develop a non-linear multivariate analysis by employing a kernel method using two homologous high-mannose-type oligosaccharides composed of ten and eleven residues as model molecules. These oligosaccharides' conformers derived from simulations were mapped into reproductive kernel Hilbert space with a positive definite function in which all required non-redundant variables for describing the oligosaccharide conformations can be treated in a non-biased manner. By applying Gaussian mixture model clustering, the oligosaccharide conformers were successfully classified by different funnels in the free-energy landscape, enabling a systematic comparison of conformational ensembles of the homologous oligosaccharides. The results shed light on the contributions of intraresidue conformational factors such as the hydroxyl group orientation and/or ring puckering state to their global conformational dynamics. Our methodology will open opportunities to explore oligosaccharides' conformational spaces, and more generally, molecules with high degrees of motional freedom.

OpenAWSEM with Open3SPN2: A fast, flexible, and accessible framework for large-scale coarse-grained biomolecular simulations.

We present OpenAWSEM and Open3SPN2, new cross-compatible implementations of coarse-grained models for protein (AWSEM) and DNA (3SPN2) molecular dynamics simulations within the OpenMM framework. These new implementations retain the chemical accuracy and intrinsic efficiency of the original models while adding GPU acceleration and the ease of forcefield modification provided by OpenMM's Custom Forces software framework. By utilizing GPUs, we achieve around a 30-fold speedup in protein and protein-DNA simulations over the existing LAMMPS-based implementations running on a single CPU core. We showcase the benefits of OpenMM's Custom Forces framework by devising and implementing two new potentials that allow us to address important aspects of protein folding and structure prediction and by testing the ability of the combined OpenAWSEM and Open3SPN2 to model protein-DNA binding. The first potential is used to describe the changes in effective interactions that occur as a protein becomes partially buried in a membrane. We also introduced an interaction to describe proteins with multiple disulfide bonds. Using simple pairwise disulfide bonding terms results in unphysical clustering of cysteine residues, posing a problem when simulating the folding of proteins with many cysteines. We now can computationally reproduce Anfinsen's early Nobel prize winning experiments by using OpenMM's Custom Forces framework to introduce a multi-body disulfide bonding term that prevents unphysical clustering. Our protein-DNA simulations show that the binding landscape is funneled towards structures that are quite similar to those found using experiments. In summary, this paper provides a simulation tool for the molecular biophysics community that is both easy to use and sufficiently efficient to simulate large proteins and large protein-DNA systems that are central to many cellular processes. These codes should facilitate the interplay between molecular simulations and cellular studies, which have been hampered by the large mismatch between the time and length scales accessible to molecular simulations and those relevant to cell biology.

BioAFMviewer: An interactive interface for simulated AFM scanning of biomolecular structures and dynamics.

We provide a stand-alone software, the BioAFMviewer, which transforms biomolecular structures into the graphical representation corresponding to the outcome of atomic force microscopy (AFM) experiments. The AFM graphics is obtained by performing simulated scanning over the molecular structure encoded in the corresponding PDB file. A versatile molecular viewer integrates the visualization of PDB structures and control over their orientation, while synchronized simulated scanning with variable spatial resolution and tip-shape geometry produces the corresponding AFM graphics. We demonstrate the applicability of the BioAFMviewer by comparing simulated AFM graphics to high-speed AFM observations of proteins. The software can furthermore process molecular movies of conformational motions, e.g. those obtained from servers which model functional transitions within a protein, and produce the corresponding simulated AFM movie. The BioAFMviewer software provides the platform to employ the plethora of structural and dynamical data of proteins in order to help in the interpretation of biomolecular AFM experiments.

Validation of DBFOLD: An efficient algorithm for computing folding pathways of complex proteins.

Atomistic simulations can provide valuable, experimentally-verifiable insights into protein folding mechanisms, but existing ab initio simulation methods are restricted to only the smallest proteins due to severe computational speed limits. The folding of larger proteins has been studied using native-centric potential functions, but such models omit the potentially crucial role of non-native interactions. Here, we present an algorithm, entitled DBFOLD, which can predict folding pathways for a wide range of proteins while accounting for the effects of non-native contacts. In addition, DBFOLD can predict the relative rates of different transitions within a protein's folding pathway. To accomplish this, rather than directly simulating folding, our method combines equilibrium Monte-Carlo simulations, which deploy enhanced sampling, with unfolding simulations at high temperatures. We show that under certain conditions, trajectories from these two types of simulations can be jointly analyzed to compute unknown folding rates from detailed balance. This requires inferring free energies from the equilibrium simulations, and extrapolating transition rates from the unfolding simulations to lower, physiologically-reasonable temperatures at which the native state is marginally stable. As a proof of principle, we show that our method can accurately predict folding pathways and Monte-Carlo rates for the well-characterized Streptococcal protein G. We then show that our method significantly reduces the amount of computation time required to compute the folding pathways of large, misfolding-prone proteins that lie beyond the reach of existing direct simulation. Our algorithm, which is available online, can generate detailed atomistic models of protein folding mechanisms while shedding light on the role of non-native intermediates which may crucially affect organismal fitness and are frequently implicated in disease.

Learning reaction coordinates via cross-entropy minimization: Application to alanine dipeptide.

We propose a cross-entropy minimization method for finding the reaction coordinate from a large number of collective variables in complex molecular systems. This method is an extension of the likelihood maximization approach describing the committor function with a sigmoid. By design, the reaction coordinate as a function of various collective variables is optimized such that the distribution of the committor pB * values generated from molecular dynamics simulations can be described in a sigmoidal manner. We also introduce the L2-norm regularization used in the machine learning field to prevent overfitting when the number of considered collective variables is large. The current method is applied to study the isomerization of alanine dipeptide in vacuum, where 45 dihedral angles are used as candidate variables. The regularization parameter is determined by cross-validation using training and test datasets. It is demonstrated that the optimal reaction coordinate involves important dihedral angles, which are consistent with the previously reported results. Furthermore, the points with pB *∼0.5 clearly indicate a separatrix distinguishing reactant and product states on the potential of mean force using the extracted dihedral angles.

EasyAmber: A comprehensive toolbox to automate the molecular dynamics simulation of proteins.

Conformational plasticity of the functionally important regions and binding sites in protein/enzyme structures is one of the key factors affecting their function and interaction with substrates/ligands. Molecular dynamics (MD) can address the challenge of accounting for protein flexibility by predicting the time-dependent behavior of a molecular system. It has a potential of becoming a particularly important tool in protein engineering and drug discovery, but requires specialized training and skills, what impedes practical use by many investigators. We have developed the easyAmber - a comprehensive set of programs to automate the molecular dynamics routines implemented in the Amber package. The toolbox can address a wide set of tasks in computational biology struggling to account for protein flexibility. The automated workflow includes a complete set of steps from the initial "static" molecular model to the MD "production run": the full-atom model building, optimization/equilibration of the molecular system, classical/conventional and accelerated molecular dynamics simulations. The easyAmber implements advanced MD protocols, but is highly automated and easy-to-operate to attract a broad audience. The toolbox can be used on a personal desktop station equipped with a compatible gaming GPU-accelerator, as well as help to manage huge workloads on a powerful supercomputer. The software provides an opportunity to operate multiple simulations of different proteins at the same time, thus significantly increasing work efficiency. The easyAmber takes the molecular dynamics to the next level in terms of usability for complex processing of large volumes of data, thus supporting the recent trend away from inefficient "static" approaches in biology toward a deeper understanding of the dynamics in protein structures. The software is freely available for download at https://biokinet.belozersky.msu.ru/easyAmber, no login required.

A multi-modal coarse grained model of DNA flexibility mappable to the atomistic level.

We present a new coarse grained method for the simulation of duplex DNA. The algorithm uses a generalized multi-harmonic model that can represent any multi-normal distribution of helical parameters, thus avoiding caveats of current mesoscopic models for DNA simulation and representing a breakthrough in the field. The method has been parameterized from accurate parmbsc1 atomistic molecular dynamics simulations of all unique tetranucleotide sequences of DNA embedded in long duplexes and takes advantage of the correlation between helical states and backbone configurations to derive atomistic representations of DNA. The algorithm, which is implemented in a simple web interface and in a standalone package reproduces with high computational efficiency the structural landscape of long segments of DNA untreatable by atomistic molecular dynamics simulations.

Revisiting the "satisfaction of spatial restraints" approach of MODELLER for protein homology modeling.

The most frequently used approach for protein structure prediction is currently homology modeling. The 3D model building phase of this methodology is critical for obtaining an accurate and biologically useful prediction. The most widely employed tool to perform this task is MODELLER. This program implements the "modeling by satisfaction of spatial restraints" strategy and its core algorithm has not been altered significantly since the early 1990s. In this work, we have explored the idea of modifying MODELLER with two effective, yet computationally light strategies to improve its 3D modeling performance. Firstly, we have investigated how the level of accuracy in the estimation of structural variability between a target protein and its templates in the form of σ values profoundly influences 3D modeling. We show that the σ values produced by MODELLER are on average weakly correlated to the true level of structural divergence between target-template pairs and that increasing this correlation greatly improves the program's predictions, especially in multiple-template modeling. Secondly, we have inquired into how the incorporation of statistical potential terms (such as the DOPE potential) in the MODELLER's objective function impacts positively 3D modeling quality by providing a small but consistent improvement in metrics such as GDT-HA and lDDT and a large increase in stereochemical quality. Python modules to harness this second strategy are freely available at https://github.com/pymodproject/altmod. In summary, we show that there is a large room for improving MODELLER in terms of 3D modeling quality and we propose strategies that could be pursued in order to further increase its performance.

Sequence specificity despite intrinsic disorder: How a disease-associated Val/Met polymorphism rearranges tertiary interactions in a long disordered protein.

The role of electrostatic interactions and mutations that change charge states in intrinsically disordered proteins (IDPs) is well-established, but many disease-associated mutations in IDPs are charge-neutral. The Val66Met single nucleotide polymorphism (SNP) in precursor brain-derived neurotrophic factor (BDNF) is one of the earliest SNPs to be associated with neuropsychiatric disorders, and the underlying molecular mechanism is unknown. Here we report on over 250 µs of fully-atomistic, explicit solvent, temperature replica-exchange molecular dynamics (MD) simulations of the 91 residue BDNF prodomain, for both the V66 and M66 sequence. The simulations were able to correctly reproduce the location of both local and non-local secondary structure changes due to the Val66Met mutation, when compared with NMR spectroscopy. We find that the change in local structure is mediated via entropic and sequence specific effects. We developed a hierarchical sequence-based framework for analysis and conceptualization, which first identifies "blobs" of 4-15 residues representing local globular regions or linkers. We use this framework within a novel test for enrichment of higher-order (tertiary) structure in disordered proteins; the size and shape of each blob is extracted from MD simulation of the real protein (RP), and used to parameterize a self-avoiding heterogenous polymer (SAHP). The SAHP version of the BDNF prodomain suggested a protein segmented into three regions, with a central long, highly disordered polyampholyte linker separating two globular regions. This effective segmentation was also observed in full simulations of the RP, but the Val66Met substitution significantly increased interactions across the linker, as well as the number of participating residues. The Val66Met substitution replaces ß-bridging between V66 and V94 (on either side of the linker) with specific side-chain interactions between M66 and M95. The protein backbone in the vicinity of M95 is then free to form ß-bridges with residues 31-41 near the N-terminus, which condenses the protein. A significant role for Met/Met interactions is consistent with previously-observed non-local effects of the Val66Met SNP, as well as established interactions between the Met66 sequence and a Met-rich receptor that initiates neuronal growth cone retraction.

ReaDDy 2: Fast and flexible software framework for interacting-particle reaction dynamics.

Interacting-particle reaction dynamics (iPRD) combines the simulation of dynamical trajectories of interacting particles as in molecular dynamics (MD) simulations with reaction kinetics, in which particles appear, disappear, or change their type and interactions based on a set of reaction rules. This combination facilitates the simulation of reaction kinetics in crowded environments, involving complex molecular geometries such as polymers, and employing complex reaction mechanisms such as breaking and fusion of polymers. iPRD simulations are ideal to simulate the detailed spatiotemporal reaction mechanism in complex and dense environments, such as in signalling processes at cellular membranes, or in nano- to microscale chemical reactors. Here we introduce the iPRD software ReaDDy 2, which provides a Python interface in which the simulation environment, particle interactions and reaction rules can be conveniently defined and the simulation can be run, stored and analyzed. A C++ interface is available to enable deeper and more flexible interactions with the framework. The main computational work of ReaDDy 2 is done in hardware-specific simulation kernels. While the version introduced here provides single- and multi-threading CPU kernels, the architecture is ready to implement GPU and multi-node kernels. We demonstrate the efficiency and validity of ReaDDy 2 using several benchmark examples. ReaDDy 2 is available at the https://readdy.github.io/ website.

A comprehensive ensemble model for comparing the allosteric effect of ordered and disordered proteins.

Intrinsically disordered proteins/regions (IDPs/IDRs) are prevalent in allosteric regulation. It was previously thought that intrinsic disorder is favorable for maximizing the allosteric coupling. Here, we propose a comprehensive ensemble model to compare the roles of both order-order transition and disorder-order transition in allosteric effect. It is revealed that the MWC pathway (order-order transition) has a higher probability than the EAM pathway (disorder-order transition) in allostery, suggesting a complicated role of IDPs/IDRs in regulatory proteins. In addition, an analytic formula for the maximal allosteric coupling response is obtained, which shows that too stable or too unstable state is unfavorable to endow allostery, and is thus helpful for rational design of allosteric drugs.

Trajectory-based training enables protein simulations with accurate folding and Boltzmann ensembles in cpu-hours.

An ongoing challenge in protein chemistry is to identify the underlying interaction energies that capture protein dynamics. The traditional trade-off in biomolecular simulation between accuracy and computational efficiency is predicated on the assumption that detailed force fields are typically well-parameterized, obtaining a significant fraction of possible accuracy. We re-examine this trade-off in the more realistic regime in which parameterization is a greater source of error than the level of detail in the force field. To address parameterization of coarse-grained force fields, we use the contrastive divergence technique from machine learning to train from simulations of 450 proteins. In our procedure, the computational efficiency of the model enables high accuracy through the precise tuning of the Boltzmann ensemble. This method is applied to our recently developed Upside model, where the free energy for side chains is rapidly calculated at every time-step, allowing for a smooth energy landscape without steric rattling of the side chains. After this contrastive divergence training, the model is able to de novo fold proteins up to 100 residues on a single core in days. This improved Upside model provides a starting point both for investigation of folding dynamics and as an inexpensive Bayesian prior for protein physics that can be integrated with additional experimental or bioinformatic data.

Accurate calculation of side chain packing and free energy with applications to protein molecular dynamics.

To address the large gap between time scales that can be easily reached by molecular simulations and those required to understand protein dynamics, we present a rapid self-consistent approximation of the side chain free energy at every integration step. In analogy with the adiabatic Born-Oppenheimer approximation for electronic structure, the protein backbone dynamics are simulated as preceding according to the dictates of the free energy of an instantaneously-equilibrated side chain potential. The side chain free energy is computed on the fly, allowing the protein backbone dynamics to traverse a greatly smoothed energetic landscape. This computation results in extremely rapid equilibration and sampling of the Boltzmann distribution. Our method, termed Upside, employs a reduced model involving the three backbone atoms, along with the carbonyl oxygen and amide proton, and a single (oriented) side chain bead having multiple locations reflecting the conformational diversity of the side chain's rotameric states. We also introduce a novel, maximum-likelihood method to parameterize the side chain interactions using protein structures. We demonstrate state-of-the-art accuracy for predicting χ1 rotamer states while consuming only milliseconds of CPU time. Our method enables rapidly equilibrating coarse-grained simulations that can nonetheless contain significant molecular detail. We also show that the resulting free energies of the side chains are sufficiently accurate for de novo folding of some proteins.

Kirkwood-Buff-Derived Alcohol Parameters for Aqueous Carbohydrates and Their Application to Preferential Interaction Coefficient Calculations of Proteins.

The CHARMM36 carbohydrate parameter set did not adequately reproduce experimental thermodynamic data of carbohydrate interactions with water or proteins or carbohydrate self-association; thus, a new nonbonded parameter set for carbohydrates was developed. The parameters were developed to reproduce experimental Kirkwood-Buff integral values, defined by the Kirkwood-Buff theory of solutions, and applied to simulations of glycerol, sorbitol, glucose, sucrose, and trehalose. Compared to the CHARMM36 carbohydrate parameters, these new Kirkwood-Buff-based parameters reproduced accurately carbohydrate self-association and the trend of activity coefficient derivative changes with concentration. When using these parameters, preferential interaction coefficients calculated from simulations of these carbohydrates and the proteins lysozyme, bovine serum albumin, α-chymotrypsinogen A, and RNase A agreed well with the experimental data, whereas use of the CHARMM36 parameters indicated preferential inclusion of carbohydrates, in disagreement with the experiment. Thus, calculating preferential interaction coefficients from simulations requires using a force field that accurately reproduces trends in the thermodynamic properties of binary excipient-water solutions, and in particular the trend in the activity coefficient derivative. Finally, the carbohydrate-protein simulations using the new parameters indicated that the carbohydrate size was a major factor in the distribution of different carbohydrates around a protein surface.

Novel Techniques for Observing Structural Dynamics of Photoresponsive Liquid Crystals.

We discuss in this article the experimental measurements of the molecules in liquid crystal (LC) phase using the time-resolved infrared (IR) vibrational spectroscopy and time-resolved electron diffraction. Liquid crystal phase is an important state of matter that exists between the solid and liquid phases and it is common in natural systems as well as in organic electronics. Liquid crystals are orientationally ordered but loosely packed, and therefore, the internal conformations and alignments of the molecular components of LCs can be modified by external stimuli. Although advanced time-resolved diffraction techniques have revealed picosecond-scale molecular dynamics of single crystals and polycrystals, direct observations of packing structures and ultrafast dynamics of soft materials have been hampered by blurry diffraction patterns. Here, we report time-resolved IR vibrational spectroscopy and electron diffractometry to acquire ultrafast snapshots of a columnar LC material bearing a photoactive core moiety. Differential-detection analyses of the combination of time-resolved IR vibrational spectroscopy and electron diffraction are powerful tools for characterizing structures and photoinduced dynamics of soft materials.