A Global Parallel Model Based Design of Experiments Method to Minimize Model Output Uncertainty

Model-based experiment design specifies the data to be collected that will most effectively characterize the biological system under study. Existing model-based design of experiment algorithms have primarily relied on Fisher Information Matrix-based methods to choose the best experiment in a sequential manner. However, these are largely local methods that require an initial estimate of the parameter values, which are often highly uncertain, particularly when data is limited. In this paper, we provide an approach to specify an informative sequence of multiple design points (parallel design) that will constrain the dynamical uncertainty of the biological system responses to within experimentally detectable limits as specified by the estimated experimental noise. The method is based upon computationally efficient sparse grids and requires only a bounded uncertain parameter space; it does not rely upon initial parameter estimates. The design sequence emerges through the use of scenario trees with experimental design points chosen to minimize the uncertainty in the predicted dynamics of the measurable responses of the system. The algorithm was illustrated herein using a T cell activation model for three problems that ranged in dimension from 2D to 19D. The results demonstrate that it is possible to extract useful information from a mathematical model where traditional model-based design of experiments approaches most certainly fail. The experiments designed via this method fully constrain the model output dynamics to within experimentally resolvable limits. The method is effective for highly uncertain biological systems characterized by deterministic mathematical models with limited data sets. Also, it is highly modular and can be modified to include a variety of methodologies such as input design and model discrimination.     

Exploring zipping and assembly as a protein folding principle

It has been proposed that proteins fold by a process called "Zipping and Assembly" (Z&A). Zipping refers to the growth of local substructures within the chain, and assembly refers to the coming together of already-formed pieces. Our interest here is in whether Z&A is a general method that can fold most of sequence space, to global minima, efficiently. Using the HP model, we can address this question by enumerating full conformation and sequence spaces. We find that Z&A reaches the global energy minimum native states, even though it searches only a very small fraction of conformational space, for most sequences in the full sequence space. We find that Z&A, a mechanism-based search, is more efficient in our tests than the replica exchange search method. Folding efficiency is increased for chains having: (a) small loop-closure steps, consistent with observations by Plaxco et al. 1998;277;985-994 that folding rates correlate with contact order, (b) neither too few nor too many nucleation sites per chain, and (c) assembly steps that do not occur too early in the folding process. We find that the efficiency increases with chain length, although our range of chain lengths is limited. We believe these insights may be useful for developing faster protein conformational search algorithms.     

Supervised cluster analysis for microarray data based on multivariate Gaussian mixture

MOTIVATION: Grouping genes having similar expression patterns is called gene clustering, which has been proved to be a useful tool for extracting underlying biological information of gene expression data. Many clustering procedures have shown success in microarray gene clustering; most of them belong to the family of heuristic clustering algorithms. Model-based algorithms are alternative clustering algorithms, which are based on the assumption that the whole set of microarray data is a finite mixture of a certain type of distributions with different parameters. Application of the model-based algorithms to unsupervised clustering has been reported. Here, for the first time, we demonstrated the use of the model-based algorithm in supervised clustering of microarray data. RESULTS: We applied the proposed methods to real gene expression data and simulated data. We showed that the supervised model-based algorithm is superior over the unsupervised method and the support vector machines (SVM) method. AVAILABILITY: The program written in the SAS language implementing methods I-III in this report is available upon request. The software of SVMs is available in the website http://svm.sdsc.edu/cgi-bin/nph-SVMsubmit.cgi     

Progress in protein-protein docking: atomic resolution predictions in the CAPRI experiment using RosettaDock with an improved treatment of side-chain flexibility

RosettaDock uses real-space Monte Carlo minimization (MCM) on both rigid-body and side-chain degrees of freedom to identify the lowest free energy docked arrangement of 2 protein structures. An improved version of the method that uses gradient-based minimization for off-rotamer side-chain optimization and includes information from unbound structures was used to create predictions for Rounds 4 and 5 of CAPRI. First, large numbers of independent MCM trajectories were carried out and the lowest free energy docked configurations identified. Second, new trajectories were started from these lowest energy structures to thoroughly sample the surrounding conformation space, and the lowest energy configurations were submitted as predictions. For all cases in which there were no significant backbone conformational changes, a small number of very low-energy configurations were identified in the first, global search and subsequently found to be close to the center of the basin of attraction in the free energy landscape in the second, local search. Following the release of the experimental coordinates, it was found that the centers of these free energy minima were remarkably close to the native structures in not only the rigid-body orientation but also the detailed conformations of the side-chains. Out of 8 targets, the lowest energy models had interface root-mean-square deviations (RMSDs) less than 1.1 A from the correct structures for 6 targets, and interface RMSDs less than 0.4 A for 3 targets. The predictions were top submissions to CAPRI for Targets 11, 12, 14, 15, and 19. The close correspondence of the lowest free energy structures found in our searches to the experimental structures suggests that our free energy function is a reasonable representation of the physical chemistry, and that the real space search with full side-chain flexibility to some extent solves the protein-protein docking problem in the absence of significant backbone conformational changes. On the other hand, the approach fails when there are significant backbone conformational changes as the steric complementarity of the 2 proteins cannot be modeled without incorporating backbone flexibility, and this is the major goal of our current work.     

Model-Based Reasoning in Humans Becomes Automatic with Training

Model-based and model-free reinforcement learning (RL) have been suggested as algorithmic realizations of goal-directed and habitual action strategies. Model-based RL is more flexible than model-free but requires sophisticated calculations using a learnt model of the world. This has led model-based RL to be identified with slow, deliberative processing, and model-free RL with fast, automatic processing. In support of this distinction, it has recently been shown that model-based reasoning is impaired by placing subjects under cognitive load—a hallmark of non-automaticity. Here, using the same task, we show that cognitive load does not impair model-based reasoning if subjects receive prior training on the task. This finding is replicated across two studies and a variety of analysis methods. Thus, task familiarity permits use of model-based reasoning in parallel with other cognitive demands. The ability to deploy model-based reasoning in an automatic, parallelizable fashion has widespread theoretical implications, particularly for the learning and execution of complex behaviors. It also suggests a range of important failure modes in psychiatric disorders.     

The Equivalence of Information-Theoretic and Likelihood-Based Methods for Neural Dimensionality Reduction

Stimulus dimensionality-reduction methods in neuroscience seek to identify a low-dimensional space of stimulus features that affect a neuron’s probability of spiking. One popular method, known as maximally informative dimensions (MID), uses an information-theoretic quantity known as “single-spike information” to identify this space. Here we examine MID from a model-based perspective. We show that MID is a maximum-likelihood estimator for the parameters of a linear-nonlinear-Poisson (LNP) model, and that the empirical single-spike information corresponds to the normalized log-likelihood under a Poisson model. This equivalence implies that MID does not necessarily find maximally informative stimulus dimensions when spiking is not well described as Poisson. We provide several examples to illustrate this shortcoming, and derive a lower bound on the information lost when spiking is Bernoulli in discrete time bins. To overcome this limitation, we introduce model-based dimensionality reduction methods for neurons with non-Poisson firing statistics, and show that they can be framed equivalently in likelihood-based or information-theoretic terms. Finally, we show how to overcome practical limitations on the number of stimulus dimensions that MID can estimate by constraining the form of the non-parametric nonlinearity in an LNP model. We illustrate these methods with simulations and data from primate visual cortex.     

Efficient Sampling in Fragment-Based Protein Structure Prediction Using an Estimation of Distribution Algorithm

Fragment assembly is a powerful method of protein structure prediction that builds protein models from a pool of candidate fragments taken from known structures. Stochastic sampling is subsequently used to refine the models. The structures are first represented as coarse-grained models and then as all-atom models for computational efficiency. Many models have to be generated independently due to the stochastic nature of the sampling methods used to search for the global minimum in a complex energy landscape. In this paper we present , a fragment-based approach which shares information between the generated models and steers the search towards native-like regions. A distribution over fragments is estimated from a pool of low energy all-atom models. This iteratively-refined distribution is used to guide the selection of fragments during the building of models for subsequent rounds of structure prediction. The use of an estimation of distribution algorithm enabled to reach lower energy levels and to generate a higher percentage of near-native models. uses an all-atom energy function and produces models with atomic resolution. We observed an improvement in energy-driven blind selection of models on a benchmark of in comparison with the AbInitioRelax protocol.     

Conformational analysis using a truth maintenance system

This paper describes Rotamer, a program we wrote to assist chemists in identifying the likely three-dimensional (3D) shapes (conformations) that a flexible molecule might assume. The search of the space of all possible conformations is approached as a constraint-satisfaction problem. The structure generator in Rotamer uses a simplified truth maintenance system (TMS) to cache information acquired in the course of exploring the conformation space. The computational benefits of using this TMS increase with the size of the problem.     

High-accuracy refinement using Rosetta in CASP13

Because proteins generally fold to their lowest free energy states, energy-guided refinement in principle should be able to systematically improve the quality of protein structure models generated using homologous structure or co-evolution derived information. However, because of the high dimensionality of the search space, there are far more ways to degrade the quality of a near native model than to improve it, and hence, refinement methods are very sensitive to energy function errors. In the 13th Critial Assessment of techniques for protein Structure Prediction (CASP13), we sought to carry out a thorough search for low energy states in the neighborhood of a starting model using restraints to avoid straying too far. The approach was reasonably successful in improving both regions largely incorrect in the starting models as well as core regions that started out closer to the correct structure. Models with GDT-HA over 70 were obtained for five targets and for one of those, an accuracy of 0.5 å backbone root-mean-square deviation (RMSD) was achieved. An important current challenge is to improve performance in refining oligomers and larger proteins, for which the search problem remains extremely difficult.     

MOLS sampling and its applications in structural biophysics

This review describes the MOLS method and its applications. This computational method has been developed in our laboratory primarily to explore the conformational space of small peptides and identify features of interest, particularly the minima, i.e., the low energy conformations. A systematic “brute-force” search through the vast conformational space for such features faces the insurmountable problem of combinatorial explosion, whilst other techniques, e.g., Monte Carlo searches, are somewhat limited in their region of exploration and may be considered inexhaustive. The MOLS method, on the other hand, uses a sampling technique commonly employed in experimental design theory to identify a small sample of the conformational space that nevertheless retains information about the entire space. The information is extracted using a technique that is a variant of the self-consistent mean field technique, which has been used to identify, for example, the optimal set of side-chain conformations in a protein. Applications of the MOLS method to understand peptide structure, predict the structures of loops in proteins, predict three-dimensional structures of small proteins, and arrive at the best conformation, orientation, and positions of a small molecule ligand in a protein receptor site have all yielded satisfactory results.     

一种在格点模型上模拟蛋白质折叠结构的优化算法

蛋白质折叠问题是生物信息学中一个经典的多项式复杂程度的非确定性(non-deterministic polynomial,NP)难度问题.势能曲面变平法(ELP)是一种启发式的全局优化算法.通过对ELP方法中的直方图函数提出一种新的更新机制,并将基于贪心策略的初始构象的产生,基于牵引移动的邻域搜索策略与ELP方法相结合,为面心立方体(FCC)格点模型的蛋白质折叠问题提出一种改进的势能曲面变平(ELP+)算法.采用文献中9条常用序列作为测试集.对于每条序列,ELP+算法均能找到与文献中的算法所得到的最低能量相等或更低的能量.实验结果表明,ELP+算法是求解FCC格点模型的蛋白质折叠问题的一种有效算法.     

Free energy landscape of protein folding in water: explicit vs. implicit solvent

The Generalized Born (GB) continuum solvent model is arguably the most widely used implicit solvent model in protein folding and protein structure prediction simulations; however, it still remains an open question on how well the model behaves in these large-scale simulations. The current study uses the beta-hairpin from C-terminus of protein G as an example to explore the folding free energy landscape with various GB models, and the results are compared to the explicit solvent simulations and experiments. All free energy landscapes are obtained from extensive conformation space sampling with a highly parallel replica exchange method. Because solvation model parameters are strongly coupled with force fields, five different force field/solvation model combinations are examined and compared in this study, namely the explicit solvent model: OPLSAA/SPC model, and the implicit solvent models: OPLSAA/SGB (Surface GB), AMBER94/GBSA (GB with Solvent Accessible Surface Area), AMBER96/GBSA, and AMBER99/GBSA. Surprisingly, we find that the free energy landscapes from implicit solvent models are quite different from that of the explicit solvent model. Except for AMBER96/GBSA, all other implicit solvent models find the lowest free energy state not the native state. All implicit solvent models show erroneous salt-bridge effects between charged residues, particularly in OPLSAA/SGB model, where the overly strong salt-bridge effect results in an overweighting of a non-native structure with one hydrophobic residue F52 expelled from the hydrophobic core in order to make better salt bridges. On the other hand, both AMBER94/GBSA and AMBER99/GBSA models turn the beta-hairpin in to an alpha-helix, and the alpha-helical content is much higher than the previously reported alpha-helices in an explicit solvent simulation with AMBER94 (AMBER94/TIP3P). Only AMBER96/GBSA shows a reasonable free energy landscape with the lowest free energy structure the native one despite an erroneous salt-bridge between D47 and K50. Detailed results on free energy contour maps, lowest free energy structures, distribution of native contacts, alpha-helical content during the folding process, NOE comparison with NMR, and temperature dependences are reported and discussed for all five models.     

Conformational sampling using high-temperature molecular dynamics

High-temperature molecular dynamics as a method for conformational search was explored on the antigen combining site of McPC 603, a phosphorylcholine binding immunoglobulin. Simulations at temperatures of 500, 800, and 1500 K were run for 111.5, 101.7, and 76.3 ps, respectively. The effectiveness of the search was assessed using a variety of methods. For the shorter hypervariable loops, molecular dynamics explored an appreciable fraction of the conformational space as evidenced by a comparison to a simple theoretical model of the size of the conformational space. However, for the longer loops and the antigen combining site as a whole, the simulation times were too short for a complete search. The simulations at 500 and 800 K both generated conformations that minimized to energies 200 kcal/mole lower than the crystal structure. However, the 1500 K simulation produced higher energy structures, even after minimization; in addition, this highest temperature run had many cis-trans peptide isomerizations. This suggests that 1500 K is too high a temperature for unconstrained conformational sampling. Comparison of the results of high temperature molecular dynamics with a direct conformational search method, [R. E. Bruccoleri & M. Karplus (1987) Biopolymers 26, 137-168]. showed that the two methods did not overlap much in conformational space. Simple geometric measures of the conformational space indicated that the direct method covered more space than molecular dynamics at the lower temperature, but not at 1500 K. The results suggest that high-temperature molecular dynamics can aid in conformational searches.     

A Fast Flexible Docking Method using an Incremental Construction Algorithm

We present an automatic method for docking organic ligands into protein binding sites. The method can be used in the design process of specific protein ligands. It combines an appropriate model of the physico-chemical properties of the docked molecules with efficient methods for sampling the conformational space of the ligand. If the ligand is flexible, it can adopt a large variety of different conformations. Each such minimum in conformational space presents a potential candidate for the conformation of the ligand in the complexed state. Our docking method samples the conformation space of the ligand on the basis of a discrete model and uses a tree-search technique for placing the ligand incrementally into the active site. For placing the first fragment of the ligand into the protein, we use hashing techniques adapted from computer vision. The incremental construction algorithm is based on a greedy strategy combined with efficient methods for overlap detection and for the search of new interactions. We present results on 19 complexes of which the binding geometry has been crystallographically determined. All considered ligands are docked in at most three minutes on a current workstation. The experimentally observed binding mode of the ligand is reproduced with 0.5 to 1.2 Å rms deviation. It is almost always found among the highest-ranking conformations computed.     

Understanding the structural characteristics of compstatin by conformational space annealing

The structural characteristics of the 13-residue compstatin molecule are investigated using the conformational space annealing (CSA) method with CHARMM force field and the GBSA continuum solvent model. In order to sample conformations in the energy range of the minimized NMR structures, we have used the stopping criterion to the CSA search when a conformation whose energy is less than -490 kcal/mol is found. With this stopping criterion, a great variety of conformations are generated around experimentally known structures. Twenty independent CSA runs starting from random states find 1000 representative conformations in the energy landscape of the compstatin, which are classified into thirty-one structural families. The majority of the conformations (94.4%) are in the coil state. Other conformers containing a 3(10)-helix, a pi-helix, a beta-hairpin, and an alpha-helix are also found.     

Conformational switching upon phosphorylation: a predictive framework based on energy landscape principles

We investigate how post-translational phosphorylation modifies the global conformation of a protein by changing its free energy landscape using two test proteins, cystatin and NtrC. We first examine the changes in a free energy landscape caused by phosphorylation using a model containing information about both structural forms. For cystatin the free energy cost is fairly large indicating a low probability of sampling the phosphorylated conformation in a perfectly funneled landscape. The predicted barrier for NtrC conformational transition is several times larger than the barrier for cystatin, indicating that the switch protein NtrC most probably follows a partial unfolding mechanism to move from one basin to the other. Principal component analysis and linear response theory show how the naturally occurring conformational changes in unmodified proteins are captured and stabilized by the change of interaction potential. We also develop a partially guided structure prediction Hamiltonian which is capable of predicting the global structure of a phosphorylated protein using only knowledge of the structure of the unphosphorylated protein or vice versa. This algorithm makes use of a generic transferable long-range residue contact potential along with details of structure short range in sequence. By comparing the results obtained with this guided transferable potential to those from the native-only, perfectly funneled Hamiltonians, we show that the transferable Hamiltonian correctly captures the nature of the global conformational changes induced by phosphorylation and can sample substantially correct structures for the modified protein with high probability.     

Explicit and implicit water simulations of a beta-hairpin peptide.

The conformational properties of a beta-hairpin peptide (YITNSDGTWT) were studied by using both explicit and implicit water simulations. The conformational space of the peptide was scanned by using a restricted hydrogen-bonding search method. The search method used generated the conformational space with enough diversity and good representation of beta-hairpin structures. By using a total surface area-based treatment of hydrophobic interactions, implicit water simulations failed to discriminate between experimental beta-hairpin structures from the rest of the conformers present in the authors' conformation library. However, with inclusion of vibrational free energy and accounting separately for polar and nonpolar surface areas, the nuclear magnetic resonance structure was ranked successfully as the most stable conformation. There is a loose correlation between the conformational energies by the continuum model and the conformational energies by explicit water simulation for conformers with similar structures. However, in terms of solvation energy, both approaches have a much better correlation. By using proper treatment of surface effect (partition of the surface area into polar and nonpolar areas) and including vibrational free-energy contribution, the continuum models should be reliable. Furthermore, the authors found that, for this peptide, beta-hairpin structures have large vibrational entropy that contributes decisively to the stability of folded beta-hairpin structures. Proteins 1999;37:73-87.     

New developments of the electrostatically driven monte carlo method: Test on the membrane-bound portion of melittin

The electrostatically driven Monte Carlo (EDMC) method has been greatly improved by adding a series of new features, including a procedure for cluster analysis of the accepted conformations. This information is used to guide the search for the global energy minimum. Alternative procedures for generating perturbed conformations to sample the conformational space were also included. These procedures enhance the efficiency of the method by generating a larger number of low-energy conformations. The improved EDMC method has been used to explore the conformational space of a 20-residue polypeptide chain whose sequence corresponds to the membrane-bound portion of melittin. The ECEPP/3 (Empirical Conformational Energy Program for Peptides) algorithm was used to describe the conformational energy of the chain. After an exhaustive search involving 14 independent runs, the lowest energy conformation (LEC) (−91.0 kcal/mol) of the entire study was encountered in four of the runs, while conformations higher in energy by no more than 1.8 kcal/mol were found in the remaining runs with the exception of one of them (run 8). The LEC is identical to the conformation found recently by J. Lee, H.A. Scheraga, and S. Rackovsky [(1998) “Conformational Analysis of the 20-Residue Membrane-Bound Portion of Melittin by Conformational Space Annealing,” Biopolymers, Vol. 46, pp. 103–115] as the lowest energy conformation obtained in their study using the conformational space annealing method. These results suggest that this conformation corresponds to the global energy minimum of the ECEPP/3 potential function for this specific sequence; it also appears to be the conformation of lowest free energy. © 1998 John Wiley & Sons, Inc. Biopoly 46: 117–126, 1998     

Iterative cluster-NMA: A tool for generating conformational transitions in proteins

Computational models provide insight into the structure-function relationship in proteins. These approaches, especially those based on normal mode analysis, can identify the accessible motion space around a given equilibrium structure. The large magnitude, collective motions identified by these methods are often well aligned with the general direction of the expected conformational transitions. However, these motions cannot realistically be extrapolated beyond the local neighborhood of the starting conformation. In this article, the iterative cluster-NMA (icNMA) method is presented for traversing the energy landscape from a starting conformation to a desired goal conformation. This is accomplished by allowing the evolving geometry of the intermediate structures to define the local accessible motion space, and thus produce an appropriate displacement. Following the derivation of the icNMA method, a set of sample simulations are performed to probe the robustness of the model. A detailed analysis of beta1,4-galactosyltransferase-T1 is also given, to highlight many of the capabilities of icNMA. Remarkably, during the transition, a helix is seen to be extended by an additional turn, emphasizing a new unknown role for secondary structures to absorb slack during transitions. The transition pathway for adenylate kinase, which has been frequently studied in the literature, is also discussed.