Integrating diffusion maps with umbrella sampling: application to alanine dipeptide

Nonlinear dimensionality reduction techniques can be applied to molecular simulation trajectories to systematically extract a small number of variables with which to parametrize the important dynamical motions of the system. For molecular systems exhibiting free energy barriers exceeding a few k(B)T, inadequate sampling of the barrier regions between stable or metastable basins can lead to a poor global characterization of the free energy landscape. We present an adaptation of a nonlinear dimensionality reduction technique known as the diffusion map that extends its applicability to biased umbrella sampling simulation trajectories in which restraining potentials are employed to drive the system into high free energy regions and improve sampling of phase space. We then propose a bootstrapped approach to iteratively discover good low-dimensional parametrizations by interleaving successive rounds of umbrella sampling and diffusion mapping, and we illustrate the technique through a study of alanine dipeptide in explicit solvent.     

Algorithmic dimensionality reduction for molecular structure analysis

Dimensionality reduction approaches have been used to exploit the redundancy in a Cartesian coordinate representation of molecular motion by producing low-dimensional representations of molecular motion. This has been used to help visualize complex energy landscapes, to extend the time scales of simulation, and to improve the efficiency of optimization. Until recently, linear approaches for dimensionality reduction have been employed. Here, we investigate the efficacy of several automated algorithms for nonlinear dimensionality reduction for representation of trans, trans-1,2,4-trifluorocyclo-octane conformation--a molecule whose structure can be described on a 2-manifold in a Cartesian coordinate phase space. We describe an efficient approach for a deterministic enumeration of ring conformations. We demonstrate a drastic improvement in dimensionality reduction with the use of nonlinear methods. We discuss the use of dimensionality reduction algorithms for estimating intrinsic dimensionality and the relationship to the Whitney embedding theorem. Additionally, we investigate the influence of the choice of high-dimensional encoding on the reduction. We show for the case studied that, in terms of reconstruction error root mean square deviation, Cartesian coordinate representations and encodings based on interatom distances provide better performance than encodings based on a dihedral angle representation.     

A normal mode-based geometric simulation approach for exploring biologically relevant conformational transitions in proteins

A three-step approach for multiscale modeling of protein conformational changes is presented that incorporates information about preferred directions of protein motions into a geometric simulation algorithm. The first two steps are based on a rigid cluster normal-mode analysis (RCNMA). Low-frequency normal modes are used in the third step (NMSim) to extend the recently introduced idea of constrained geometric simulations of diffusive motions in proteins by biasing backbone motions of the protein, whereas side-chain motions are biased toward favorable rotamer states. The generated structures are iteratively corrected regarding steric clashes and stereochemical constraint violations. The approach allows performing three simulation types: unbiased exploration of conformational space; pathway generation by a targeted simulation; and radius of gyration-guided simulation. When applied to a data set of proteins with experimentally observed conformational changes, conformational variabilities are reproduced very well for 4 out of 5 proteins that show domain motions, with correlation coefficients r > 0.70 and as high as r = 0.92 in the case of adenylate kinase. In 7 out of 8 cases, NMSim simulations starting from unbound structures are able to sample conformations that are similar (root-mean-square deviation = 1.0-3.1 ?) to ligand bound conformations. An NMSim generated pathway of conformational change of adenylate kinase correctly describes the sequence of domain closing. The NMSim approach is a computationally efficient alternative to molecular dynamics simulations for conformational sampling of proteins. The generated conformations and pathways of conformational transitions can serve as input to docking approaches or as starting points for more sophisticated sampling techniques.     

CAMERRA: An analysis tool for the computation of conformational dynamics by evaluating residue–residue associations

A computational method which extracts the dominant motions from an ensemble of biomolecular conformations via a correlation analysis of residue–residue contacts is presented. The algorithm first renders the structural information into contact matrices, then constructs the collective modes based on the correlated dynamics of a selected set of dynamic contacts. Associated programs can bridge the results for further visualization using graphics software. The aim of this method is to provide an analysis of conformations of biopolymers from the contact viewpoint. It may assist a systematical uncovering of conformational switching mechanisms existing in proteins and biopolymer systems in general by statistical analysis of simulation snapshots. In contrast to conventional correlation analyses of Cartesian coordinates (such as distance covariance analysis and Cartesian principal component analysis), this program also provides an alternative way to locate essential collective motions in general. Herein, we detail the algorithm in a stepwise manner and comment on the importance of the method as applied to decoding allosteric mechanisms. © 2018 Wiley Periodicals, Inc.     

Nonlinear mapping networks

Among the many dimensionality reduction techniques that have appeared in the statistical literature, multidimensional scaling and nonlinear mapping are unique for their conceptual simplicity and ability to reproduce the topology and structure of the data space in a faithful and unbiased manner. However, a major shortcoming of these methods is their quadratic dependence on the number of objects scaled, which imposes severe limitations on the size of data sets that can be effectively manipulated. Here we describe a novel approach that combines conventional nonlinear mapping techniques with feed-forward neural networks, and allows the processing of data sets orders of magnitude larger than those accessible with conventional methodologies. Rooted on the principle of probability sampling, the method employs a classical algorithm to project a small random sample, and then "learns" the underlying nonlinear transform using a multilayer neural network trained with the back-propagation algorithm. Once trained, the neural network can be used in a feed-forward manner to project the remaining members of the population as well as new, unseen samples with minimal distortion. Using examples from the fields of image processing and combinatorial chemistry, we demonstrate that this method can generate projections that are virtually indistinguishable from those derived by conventional approaches. The ability to encode the nonlinear transform in the form of a neural network makes nonlinear mapping applicable to a wide variety of data mining applications involving very large data sets that are otherwise computationally intractable.     

Detecting coupled collective motions in protein by independent subspace analysis

Protein dynamics evolves in a high-dimensional space, comprising aharmonic, strongly correlated motional modes. Such correlation often plays an important role in analyzing protein function. In order to identify significantly correlated collective motions, here we employ independent subspace analysis based on the subspace joint approximate diagonalization of eigenmatrices algorithm for the analysis of molecular dynamics (MD) simulation trajectories. From the 100 ns MD simulation of T4 lysozyme, we extract several independent subspaces in each of which collective modes are significantly correlated, and identify the other modes as independent. This method successfully detects the modes along which long-tailed non-Gaussian probability distributions are obtained. Based on the time cross-correlation analysis, we identified a series of events among domain motions and more localized motions in the protein, indicating the connection between the functionally relevant phenomena which have been independently revealed by experiments.     

The reweighted path ensemble

We introduce a reweighting scheme for the path ensembles in the transition interface sampling framework. The reweighting allows for the analysis of free energy landscapes and committor projections in any collective variable space. We illustrate the reweighting scheme on a two dimensional potential with a nonlinear reaction coordinate and on a more realistic simulation of the Trp-cage folding process. We suggest that the reweighted path ensemble can be used to optimize possible nonlinear reaction coordinates.     

The ellipsoid algorithm as a method for the determination of polypeptide conformations from experimental distance constraints and energy minimization

A new method for constrained nonlinear optimization known as the ellipsoid algorithm is evaluated as a means of determining and refining the conformations of peptides. Advantages of the ellipsoid algorithm over conventional optimization methods include that it avoids many local minima that other methods would be trapped by, and that it is sometimes able to find optimum solutions in which the constraints are satisfied exactly. The dihedral angles about single bonds were used as variables to keep the dimensionality low (the rate of convergence decreases rapidly with increasing dimensionality of the problem). The method is evaluated on problems involving distance constraints, and for minimization of conformational energy functions. In an initial application, conformations consistent with an experimental set of NMR distance constraints were obtained in a problem involving 48 variable dihedral angles.     

Multidimensional scaling and visualization of large molecular similarity tables

Multidimensional scaling (MDS) is a collection of statistical techniques that attempt to embed a set of patterns described by means of a dissimilarity matrix into a low‐dimensional display plane in a way that preserves their original pairwise interrelationships as closely as possible. Unfortunately, current MDS algorithms are notoriously slow, and their use is limited to small data sets. In this article, we present a family of algorithms that combine nonlinear mapping techniques with neural networks, and make possible the scaling of very large data sets that are intractable with conventional methodologies. The method employs a nonlinear mapping algorithm to project a small random sample, and then “learns” the underlying transform using one or more multilayer perceptrons. The distinct advantage of this approach is that it captures the nonlinear mapping relationship in an explicit function, and allows the scaling of additional patterns as they become available, without the need to reconstruct the entire map. A novel encoding scheme is described, allowing this methodology to be used with a wide variety of input data representations and similarity functions. The potential of the algorithm is illustrated in the analysis of two combinatorial libraries and an ensemble of molecular conformations. The method is particularly useful for extracting low‐dimensional Cartesian coordinate vectors from large binary spaces, such as those encountered in the analysis of large chemical data sets. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 488–500, 2001     

Non-Boltzmann sampling and Bennett's acceptance ratio method: how to profit from bending the rules

The exact computation of free energy differences requires adequate sampling of all relevant low energy conformations. Especially in systems with rugged energy surfaces, adequate sampling can only be achieved by biasing the exploration process, thus yielding non-Boltzmann probability distributions. To obtain correct free energy differences from such simulations, it is necessary to account for the effects of the bias in the postproduction analysis. We demonstrate that this can be accomplished quite simply with a slight modification of Bennett's Acceptance Ratio method, referring to this technique as Non-Boltzmann Bennett. We illustrate the method by several examples and show how a creative choice of the biased state(s) used during sampling can also improve the efficiency of free energy simulations.     

Free Energy Perturbation Hamiltonian Replica-Exchange Molecular Dynamics (FEP/H-REMD) for Absolute Ligand Binding Free Energy Calculations

Free Energy Perturbation with Replica Exchange Molecular Dynamics (FEP/REMD) offers a powerful strategy to improve the convergence of free energy computations. In particular, it has been shown previously that a FEP/REMD scheme allowing random moves within an extended replica ensemble of thermodynamic coupling parameters "lambda" can improve the statistical convergence in calculations of absolute binding free energy of ligands to proteins [J. Chem. Theory Comput. 2009, 5, 2583]. In the present study, FEP/REMD is extended and combined with an accelerated MD simulations method based on Hamiltonian replica-exchange MD (H-REMD) to overcome the additional problems arising from the existence of kinetically trapped conformations within the protein receptor. In the combined strategy, each system with a given thermodynamic coupling factor lambda in the extended ensemble is further coupled with a set of replicas evolving on a biased energy surface with boosting potentials used to accelerate the inter-conversion among different rotameric states of the side chains in the neighborhood of the binding site. Exchanges are allowed to occur alternatively along the axes corresponding to the thermodynamic coupling parameter lambda and the boosting potential, in an extended dual array of coupled lambda- and H-REMD simulations. The method is implemented on the basis of new extensions to the REPDSTR module of the biomolecular simulation program CHARMM. As an illustrative example, the absolute binding free energy of p-xylene to the nonpolar cavity of the L99A mutant of T4 lysozyme was calculated. The tests demonstrate that the dual lambda-REMD and H-REMD simulation scheme greatly accelerates the configurational sampling of the rotameric states of the side chains around the binding pocket, thereby improving the convergence of the FEP computations.     

On the precision of quasi steady state assumptions in stochastic dynamics

Many biochemical networks have complex multidimensional dynamics and there is a long history of methods that have been used for dimensionality reduction for such reaction networks. Usually a deterministic mass action approach is used; however, in small volumes, there are significant fluctuations from the mean which the mass action approach cannot capture. In such cases stochastic simulation methods should be used. In this paper, we evaluate the applicability of one such dimensionality reduction method, the quasi-steady state approximation (QSSA) [L. Menten and M. Michaelis, "Die kinetik der invertinwirkung," Biochem. Z 49, 333369 (1913)] for dimensionality reduction in case of stochastic dynamics. First, the applicability of QSSA approach is evaluated for a canonical system of enzyme reactions. Application of QSSA to such a reaction system in a deterministic setting leads to Michaelis-Menten reduced kinetics which can be used to derive the equilibrium concentrations of the reaction species. In the case of stochastic simulations, however, the steady state is characterized by fluctuations around the mean equilibrium concentration. Our analysis shows that a QSSA based approach for dimensionality reduction captures well the mean of the distribution as obtained from a full dimensional simulation but fails to accurately capture the distribution around that mean. Moreover, the QSSA approximation is not unique. We have then extended the analysis to a simple bistable biochemical network model proposed to account for the stability of synaptic efficacies; the substrate of learning and memory [J. E. Lisman, "A mechanism of memory storage insensitive to molecular turnover: A bistable autophosphorylating kinase," Proc. Natl. Acad. Sci. U.S.A. 82, 3055-3057 (1985)]. Our analysis shows that a QSSA based dimensionality reduction method results in errors as big as two orders of magnitude in predicting the residence times in the two stable states.     

Low-frequency localized oscillations of the DNA double strand

The results of the molecular dynamic study on the low-energy elementary linear and nonlinear excitations of DNA macromolecule obtained within the framework of the coarse-grain model of the DNA double strand are presented. The characteristics of the basic states of the model agree well with the experimental parameters of both A and B conformations of the DNA double strand. The correlation between the directly calculated dispersion curves and density distribution of the frequency spectrum obtained in the simulation experiments is found. Special attention is focused on the soliton type of nonlinear localized excitations (breathers). These excitations are shown to exist in several frequency gaps in which the propagation of harmonic linear waves is prohibited. The types of motions corresponding to all calculated breathers are identified. The correlation between the two types of breathers and their analogs studied in terms of unidimensional models and treated as elementary excitations responsible for the initial stage of the opening of the DNA double strand is established.     

Functionally relevant protein motions: extracting basin-specific collective coordinates from molecular dynamics trajectories

Functionally relevant motion of proteins has been associated with a number of atoms moving in a concerted fashion along so-called "collective coordinates." We present an approach to extract collective coordinates from conformations obtained from molecular dynamics simulations. The power of this technique for differentiating local structural fluctuations between classes of conformers obtained by clustering is illustrated by analyzing nanosecond-long trajectories for the response regulator protein Spo0F of Bacillus subtilis, generated both in vacuo and using an implicit-solvent representation. Conformational clustering is performed using automated histogram filtering of the inter-C(alpha) distances. Orthogonal (varimax) rotation of the vectors obtained by principal component analysis of these interresidue distances for the members of individual clusters is key to the interpretation of collective coordinates dominating each conformational class. The rotated loadings plots isolate significant variation in interresidue distances, and these are associated with entire mobile secondary structure elements. From this we infer concerted motions of these structural elements. For the Spo0F simulations employing an implicit-solvent representation, collective coordinates obtained in this fashion are consistent with the location of the protein's known active sites and experimentally determined mobile regions.     

Highly entangled polymer primitive chain network simulations based on dynamic tube dilation

The concept of dynamic tube dilation (DTD) is here used to formulate a new simulation scheme to obtain the linear viscoelastic response of long chains with a large number of entanglements. The new scheme is based on the primitive chain network model previously proposed by some of the authors, and successfully employed to simulate linear and nonlinear behavior of moderately entangled polymers. Scaling laws are generated by the DTD concept, and allow for prediction of the linear response of very long chains on the basis of suitable simulations performed on shorter ones, without introducing adjustable parameters. Tests of the method against existing data for linear monodisperse polyisoprene and polystyrene show good quantitative agreement.     

Supervised neighborhood preserving embedding for feature extraction and its application for soft sensor modeling

Neighborhood preserving embedding (NPE) is a useful tool for learning the manifold of high‐dimensional data. As a linear approximation of nonlinear locally linear embedding, NPE can be applied to dimensionality reduction by neighborhood preserving. However, the original NPE algorithm is an unsupervised method, which extracts features without any reference to the output information. In this paper, a supervised NPE framework is proposed for output‐related feature extraction in soft sensor applications. In the supervised NPE framework, the output information is utilized to guide the procedures for constructing the adjacent graph and calculating the weight matrix, with which the intrinsic structure of the data can be better described. For performance evaluation of the proposed method, experiments on a numerical example and an industrial debutanizer column process are carried out. The results show the effectiveness of the proposed framework.