Introduction

Molecular simulations have come to play a central role in structural biology and biophysics; they are beginning to be used in cell and systems biology as well. Simulations help us develop our intuition about the behavior of models which link biological structures to function. Protein folding, molecular recognition and ligand binding, biological machines used for transport and signaling, these are some of the research areas that have been greatly enriched by computational approaches based on molecular simulations. The design of effective potentials for modeling solvation effects implicitly, the development of new sampling methods based on replica exchange molecular dynamics, and the exploration of rare events using network models are themes in the Levy group which run through our current research highlighted below.

New multi-scale models, effective potentials, and sampling methods for molecular simulations

Implicit Solvent Models

Water plays a fundamental role in virtually all biological processes. The accurate modeling of hydration thermodynamics is therefore essential for studying protein conformational equilibria, aggregation, and binding. Explicit solvent models provide the most detailed description of hydration phenomena, but they have inherent limitations that motivate the search for other ways to represent solvation. Implicit solvent models, which are based on the statistical mechanics concept of the solvent potential of mean force, are very useful alternatives to explicit solvation for modeling protein folding and binding.

We have developed an implicit solvent effective potential (AGBNP) that is suitable for molecular dynamics simulations and high- resolution modeling. It is based on a novel implementation of the pairwise descreening Generalized Born model for the electrostatic component and a nonpolar hydration free energy estimator. The model is fully analytical with first derivatives and is computationally efficient and has been incorporated into the IMPACT molecular simulation program.

Replica Exchange Molecular Dynamics

One of the key challenges in the computer simulation of proteins at the atomic level is the sampling of conformational space. The efficiency of many common sampling protocols such as Monte-Carlo (MC) and Molecular Dynamics (MD) is limited by the need to cross high free- energy barriers and rugged energy landscapes. In the Replica Exchange (RE) algorithm many coupled simulations are run in parallel; the coupling is achieved by exchanging either thermodynamic or Hamiltonian parameters. We are exploring novel implementations of replica exchange to accelerate the convergence of protein folding simulations, and the calculation of protein-ligand binding free energies. In the latter example, the REMD coupling parameter is the protein-ligand interaction energy.

Network Models, Rare Events, and Transition Paths

We can think of the data generated by replica exchange simulations of protein folding and protein-ligand binding as generating a trace through phase space. We can organize the data using graph theoretic ideas, where the nodes are conformations and the edges represent possible jumps between closely related nodes. We use the magic of histogram re-weighting to assign relative probabilities to the jumps (edges) which connect the nodes. Furthermore, we can construct transition paths on the graphs which correspond to rare transitions between stable states, and then use the tools of transition path theory to analyze the transition path ensemble, in order to determine the number of truly different important paths and their fluxes. This is a form of multi-scale modeling where the underlying data is derived from atomic simulations, while the transition paths are constructed ex post facto.

Protein folding and misfolding: fundamental physics, health-related applications

Kinetic Network Models and Protein Folding

Protein folding is a fundamental problem in modern molecular biophysics and is an example of a slow process occurring via rare events in a high-dimensional configurational space. For this reason, it is difficult for an all-atom simulation to obtain meaningful information on the kinetics and pathways of such processes. A number of strategies for addressing this problem have been proposed over the years that involve focusing on the important slow processes while neglecting the less interesting rapid kinetics by simplification of the state space or reduction of dimensionality. Generalized ensemble methods such as replica exchange molecular dynamics (REMD) have been developed that enhance the ability to obtain accurate canonical populations in complex systems by increasing sampling efficiency. However, since REMD involves temperature swaps between MD trajectories, it is not straightforward to obtain kinetic information from such simulations. We have chosen to make use of a kinetic network model in which the nodes correspond to discrete molecular conformations from REMD simulation trajectories (rather than macrostates), and the edges are derived from an ansatz based on structural similarity. Since the network is a discretized representation of the system and does not require additional energy and force evaluations, there is a considerable gain in efficiency, allowing us to study much slower kinetic processes than would be accessible using conventional MD. Since the network topology can be constructed based on virtually all degrees of freedom, this allows for multiple pathways and complex transition states. Application of the kinetic network model to study the folding of mini-proteins such as Trp-Cage, and to the conformational reorganization of the HIV-1 Protease receptor pocket, provides us with detailed information about the rich variety of transition pathways both their structural characteristics and their kinetics.

Health-related applications: α synuclein and Parkinson's disease

Unlike globular proteins with stable native structures, partially folded proteins, and natively unfolded proteins are best characterized as conformational ensembles of rapidly interconverting structures. Structural characterization of these dynamic systems is critical to understanding the basis for protein misfolding that results in protein aggregation and disease. In collaboration with Jean Baum's research group we are integrating experiments and computational models to develop a molecular picture of the steps involved in the aggregation of a synuclein, an intrinsically disordered protein that appears in aggregated forms in the brains of patients with Parkinson's disease. The conversion from the monomeric disordered state to the highly ordered aggregate fibril form is both complex and not well understood. Aggregation rates are sensitive to changes in amino acid sequence and environmental conditions. We are focusing on the early steps in the aggregation process.

The energy landscape for ligand binding to protein receptors and drug design

Understanding how to accurately predict binding free energies of small molecule ligands to protein receptor targets is a very active area of research in our lab. The statistical mechanics of binding can be formulated in many alternative ways; we are pursuing approaches particularly well suited to modern cluster computing environments with highly parallelized computation.

Statistical mechanics and the free energy of binding

A new method has been developed for estimating protein-ligand affinities with implicit solvation called the Binding Energy Distribution Analysis Method (BEDAM) which is a statistical mechanics formulation for the free energy of binding based on the analysis of the distribution of binding energies between the ligand and the protein receptor obtained in a fictitious ensemble in which the ligand resides in the binding pocket without interacting with the receptor. This framework is developed based on implicit solvation, Hamiltonian replica exchange parallel molecular dynamics sampling, and reweighting techniques. The formulation takes into account ligand-receptor interactions, multiple binding modes, and conformational entropy losses and is currently being applied to a series of important medicinal targets.

Receptor reorganization in protein-ligand binding: HIV-1 Protease and Reverse Transcriptase

We are simultaneously investigating computational approaches to calculate receptor strain energies that are a component of binding free energies. The overall affinity of a ligand for a receptor can be expressed as a balance between the strength of the interactions of the ligand for a particular binding-competent conformation of the receptor and the probability of occurrence of that conformation in the absence of the ligand. For receptors that do not experience much conformational change upon binding, the affinity is primarily based on the interaction between the ligand and receptor. However, for flexible receptors such as HIV-1 reverse transcriptase (RT), HIV-1 protease and a number of kinases, the conformational changes in the receptor may require inclusion of receptor reorganization or strain energy to properly model the binding of ligands to that protein. We are exploring the use of bioinformatic techniques and simulations to model the receptor strain energy for flexible receptors.

We are also working with collaborators at Rutgers (Eddy Arnold group) and the University of Pittsburgh School of Medicine (Mike Parniak group) to identify lead compounds that will inhibit the ribonuclease H activity of HIV-1 RT. Structure based design of inhibitors for the RNase H function of RT has lagged behind inhibitor design for the polymerase function of RT due to limited structural information about the binding modes of inhibitors to the RNase H active site. Our current strategy attempts to leverage experimental information generated by our collaborators which has identified active compounds in large libraries by comparing this information with the predictions of enrichment for the same compounds when they are docked to many putative receptor sites on the RNase H domain.

Protein stability and the evolution of drug resistance: bioinformatics and biophysics

Mutational patterns in HIV-1 Protease

Proteins evolve through random mutagenesis and their evolutionary selection is constrained by structural, functional and environmental limitations. Thermodynamic stability is by far the most important structural factor, as most proteins need to be folded in order to function. The stability range of proteins is narrow, and is estimated experimentally to be ~ 10 kcal/mol or less. Proteins operate "on a knife's edge," whereby a single highly deleterious mutation could potentially lead to an unfolded protein. Drug resistance is acquired through mutations. Primary mutations are typically destabilizing, and must be accompanied by or correlated with compensatory stabilizing mutations. We are developing statistical methods to model observed mutational patterns within proteins evolving in response to drug pressure. The system we are currently working on is HIV-1 protease for which we have access to over 40,000 publicly available amino acid sequences, isolated from patients undergoing chemotherapy. The mutational patterns in HIV protease are highly non-trivial and involve primary mutations that confer drug resistance, compensatory stabilizing effects, and viral evolvability. Our aim is to rationalize the mutational patterns using coarse grained biophysical models and to thereby relate them to the structure and energetics of the enzyme.

Computational mutagenesis methods, ranging from statistical, empirical, and physical approaches, have been useful for understanding and predicting protein stabilities. Most existing methods are limited because they are unable to reproduce correctly the magnitude of the free energy change (ddG), the difference in free energy of unfolding between wild-type and mutant proteins (> 1.0 kcal/mol deviation from experimental ddG values). In order to overcome this problem, we are developing a computational approach that uses sampling combined with a linear interaction energy (LIE) method to predict the changes in the free energy of the native state induced by mutations. Initial tests of the method have shown an unsigned error between calculated and experimental values of <1 kcal/mol. Further development of the LIE functional form in combination with sampling techniques is an active, promising project which will allow for more accurate prediction of mutational effects on protein stability which are important to the understanding of drug-resistant mutations that are formed in many protein targets including HIV-1 protease.