Robert I. Cukier

Research Interests

We use the methods of Statistical Mechanics and Quantum Mechanics to create theories of and carry out simulations on:


Markov Programs


·       README.TXT

·       analyticEns.cpp

·       DihedToCC.cpp

·       nr3.h

·       svd.h

Methods Development

CUKMODY: We have developed CUKMODY, a molecular dynamics code oriented toward simulating proteins with explicit solvation. CUKMODY (CUKierMOlecularDYnamics) is written in C++, because of issues of code reusability, modularity, readability, and maintainability. It permits the inclusion of new molecule types, is explicitly designed for large systems, uses the particle mesh Ewald (PME) method to deal with long-range electrostatic interactions, and is parallelized using pthreads so it runs on any platform. Accelerating the exploration of configuration space is an outstanding problem in MD simulation. Our methods discussed below for overcoming barriers in the high dimensional configuration space of proteins are fully integrated into CUKMODY.

The coupling of MD to quantum mechanics for our studies on proton transfer/transport is also implemented in CUKMODY.

MD is plagued by the sampling problem. Barriers in configuration space necessitate long, expensive trajectories. We develop methods that can accelerate the exploration of configuration space and permit sampling of various functional motions. These methods are variants of the Replica Exchange Method (REM) as well as various reweighting techniques. They all have the feature of speeding up the transitions over energy barriers in the high dimension potential energy describing proteins. When a reaction coordinate is known then distance REM (DREM) is useful. When parts of the Hamiltonian are modified to enchance sampling then Hamiltonian REM (HREM) can be used. And a number of variants and combinations of DREM and HREM are useful. All these methods are implementd in CUKMODY.

ANALYZER: With former graduate student Dr. H. Lou, we developed ANALYZER, a versatile trajectory analysis program that, in addition to RMSF, RMSD, H-bond, electrostatic and salt-bridge analysis, has numerous correlation and Principal Component Analysis methods. Written in C++ to take advantage of tools such as the Standard Template Library, it has a lexical parser that makes it simple to select (not select) atoms/types of atoms/all atoms of a class/. All functionality is described via "help topic" and we have had several groups at MSU use it with satisfaction after minimal training.

CCA/MCA: Trajectory analysis of proteins and other complex materials is difficult because of the desire to isolate important, functional motions from the background "noise". Principal Component Analysis is one by now standard method for extracting functional motions from trajectory data. We are extending other statistics based methods such as Canonical Correlation Analysis (CCA) and Maximum Covariance Analysis (MCA) that have as their goals ferreting out the potential correlations between different sets of variables in a trajectory. For example, the potential for correlations between the backbone and side chain atoms in a protein. We are implementing these methods in the context of protein trajectory data that require additional analytic developments.

STATE COUNTING: We have developed C++ programs for state counting in proteins. We use the dihedral rotamers of a protein backbone (the phi and psi angles) to provide a discrete basis of protein conformations. Using this discrete basis we analyze well-sampled MD trajectories to obtain the states and their populations and the corresponding conformational entropy. The state space is exponentially large and useful aspects of C++ are employed to be able to deal with state spaces of the order of 100,000,000 conformations.


Biological function relies on the chemical reactions of electron and proton transfer that take place in enzyme active centers and are strongly influenced by the geometry and energetics of the surrounding medium. The small masses of electrons and protons necessitates a quantum mechanical treatment and, as electrons and protons are charged species, they strongly interact with their surroundings, the many strongly polar amino acids in the surrounding proteins. It is now appreciated that a coupling of proton motion to electron transfer is basic to mechanisms of biological energy conversion. We develop theories to predict rates of such reactions. For PCET, we focus on how to formulate rate constant expressions for the concerted transfer of an electron and a proton and how the rate of this process compares with that of the consecutive transfer of an electron and a proton [PDF]. Stimulated by these theories, we develop computer programs to carry out quantum molecular dynamics simulations that couple the quantum electron and/or proton to the surrounding classically treated medium [PDF].

Chains of hydrogen bonded water and/or amino acid residues are essential to proton translocation, the movement of protons across membranes to create electrochemical gradients for energy transduction. We have developed methods to treat the protons involved in translocation quantum mechanically and have coupled them to the surrounding protein in order to simulate rates of this process [PDF]. The methods have been applied to cytochrome c oxidase [PDF], an enzyme responsible for converting the energy from redox and oxygen chemistry to support proton translocation and create a transmembrane pH gradient used for ATP production. Important to oxidase function is the formation of hydrogen bonded chains of waters, which can translocate protons to various active centers in the enzyme. Methods have been developed to construct and identify such chains within proteins and used to identify them in molecular dynamics simulations [PDF].

Protein stability is based on a delicate balance between energetic and entropic factors. Intrinsically disordered proteins (IDPs) interacting with a folded partner protein in the act of binding can order the IDP to form the correct functional interface by decrease in the overall free energy. We evaluate the part of the entropic cost of ordering an IDP arising from their dihedral states.  The first IDP studied [PDF] is a leucine zipper dimer that simulated with Molecular Dynamics. It does show disorder in six phi and psi dihedral angles of the N terminal sequence of one monomer. Clustering is carried out with a version of a k-means algorithm that accounts for the circular nature of dihedral angles. For the twelve dihedrals each found to have three conformers, among the resulting 531,441 states, their populations show that the first 100 (500) most populated states account for ~65% (~90%)  of the entire population, indicating that there are strong dependencies among the dihedrals’ conformations. These state populations are used to evaluate a Kullback-Leibler divergence entropy measure and obtain the dihedral configurational entropy S.  At 300 K, TS~3 kcal/mol, showing that IDP entropy, while roughly half that would be expected from independently distributed dihedrals, can be a decisive contributor to the free energy of this IDP binding and ordering

Maintaining high fidelity during transcription is essential for the accurate transfer of genetic information from DNA to RNA. RNA polymerase (RNAP) misincorporates only one wrong nucleotide per100 000 bases. We study RNAP using MD[PDF] . It is an enormous molecular machine with many residues, DNA, RNA and NTPs. Using MD simulations, the effects of different active site NTPs in both open and closed trigger loop (TL) structures of RNAPs are compared. Unfavorable initial binding of mismatched substrates in the active site with an open TL is proposed to be the first fidelity checkpoint. The leaving of an incorrect substrate is much easier than a correct one energetically from the umbrella sampling simulations. Then, the closing motion of the TL, required for catalysis, is hindered by the presence of mismatched NTPs. Mismatched NTPs also lead to conformational changes in the active site, which perturb the coordination of magnesium ions and likely affect the ability to proceed with catalysis. This step appears to be the most important checkpoint for deoxy-NTP discrimination. Finally, structural perturbations in the template DNA and the nascent RNA in the presence of mismatches likely hinder nucleotide addition and provide the structural foundation for backtracking followed by removing erroneously incorporated nucleotides during proofreading.

Another promising method of analysis of MD based data is via the construction of Markov State Models (MSMs). In this approach, the MD trajectory data is clustered into conformations based on some clustering criterion – dihedral angles for example are very useful for peptides and proteins. Once states have been identified issues such as the pathways between some initial and final conformation through a number of intermediates can be identified using methods such as transition path theory (TPT). A first study[PDF] constructed a MSM for a five residue peptide, met-enkephalin, where least costly pathways between a closed and open form of the peptide are found. Even in this small system there is a rich structure of intermediate states and reasons for why some pathways are favored can be identified. Another study[PDF] on a response regulator (RR) involved in  biological signal transduction involved extensive targeted MD simulations to span inactive and active forms of a RR. MSM construction  provides a set of intermediates and TPT analysis suggests important pathways between active and inactive forms via intermediates with functional significance.  

Many proteins undergo large scale conformational changes when they bind ligands responsible for their catalytic activity. For example, Adenylate Kinase (AK) will bind ATP and AMP and, as it does so, the protein closes over these substrates to form a catalytically active complex, with water excluded. The reaction is phosphate transfer to form two ADPs, and we have studied the ternary complex formed from ATP, AMP and AK [PDF]. Exploring the closing motion of proteins with conventional Molecular Dynamics is unrealistic on normal simulation time scales. We have investigated what information from normal MD simulations at ambient and elevated temperatures can be obtained with the use of Principal Component Analysis (PCA), a method that resolves trajectory data into its important fluctuational components [PDF]. We develop methods that can accelerate the exploration of configuration space and permit sampling of the motion along the path from open to closed-like protein configurations [PDF]. These methods are variants of the Replica Exchange Method as well as various reweighting techniques [PDF]. They all have the feature of speeding up the transitions over energy barriers in the high dimension potential energy describing proteins. Applied to Adenylate Kinase [PDF], the potential of mean force along a closing reaction coordinate (RC) of about 30 A extent was successfully simulated and compared with FRET-based experiments. Remarkably, the potential of mean force has a distinct minimum at an RC value that is close to the ligand bound geometry. Therefore, apo AK can readily fluctuate to conformations favorable for ligand binding. HIV protease is a dimeric protein that also undergoes significant conformational changes in binding its ligands. Long MD simulations of a wild type (WT) and multi drug resistant (MDR) variety were carried out and analyzed with the help of PCA to contrast the motions of the WT and MDR versions of HIV protease [PDF].

Enzymes function by assembling exquisitly designed active sites with their substrates properly aligned and interacting with key residues that are responsible for the Chemistry. To establish mechanisms, the active site should be treated using Quantum Chemical methods, and these need to be at a sufficiently high level to produce reliable results. Often, the active site Chemistry is strongly influenced by the surrounding residues and these effects can be accounted for with lower-level Quantum Chemical methods and/or Molecular Mechanical force fields. We use the ONIOM method introduced by Morokuma for the former QM/QM studies and have developed QM/MM methods for the latter approach. The reaction path for the deamination reaction catalyzed by yeast cytosine deaminase (yCD), a zinc metalloenzyme of significant biomedical interest, was investigated using the ONIOM method [PDF], as well as a combined ONIOM and Molecular Dynamics method [PDF]. Guanine deaminase was also studied by similar methods [PDF].

The properties and modes of recognition of physiological DNAs associated with the four natural nucleobases can be extended, in principle, by the design of non-natural nucleobase derivatives. The goal is expansion of the genetic alphabet, with the possible outcome of producing new DNAs that have improved physical or biological properties. In collaboration with Professor Y. Bu (Institute of Theoretical Chemistry, Shandong University, PRC), hetero-ring-expanded base analogs are designed, and their structural characteristics and electronic properties are determined by density functional theory and their stability determined by molecular dynamics. The simulations we carry out show that the designed motifs can form stable DNA-like structures with improved electronic properties. For example, the HOMO-LUMO energy gaps for most of the size-expanded guanine analogs and their Watson Crick base pairs are considerably lower than those of the corresponding natural base and base pairs [PDF].

The nature of an excess electron (EE) in various media has been extensively explored because of its fundamental importance and relevance to a large class of physical and chemical phenomena associated with charge migration and radical reactions. EEs have been shown to exist in either solvated and/or surface states in some molecular solvents such as water and alcohols, and in a localized F-center-like state in alkali halide molten salts. Recently developed room-temperature ionic liquids (IL) appear to be promising “green” media and have found many intriguing applications in synthetic chemistry, separation science, materials, and electrochemical devices such as batteries and solar cells. EEs in ILs offer the possibility of new and improved reactivity. In collaboration with Professor Y. Bu (Institute of Theoretical Chemistry, Shandong University, PRC), we are using ab initio molecular dynamics (AIMD) methods to explore the states and evolution dynamics of migrating electrons in ILs [PDF].


Center for Biological Modeling

I am a founding member and former director of the Center for Biological Modeling (CBM).

The CBM (now the QBI - Quantitative Biology Initiative) brings together faculty and students from the physical, biological and computer sciences to work on problems of biological interest. There is support for graduate students interested in pursuing interdisciplinary research in this area.

Some quotes I like







Department of Chemistry
Michigan State University
East Lansing, MI 48824-1322


(517) 353-1170


( )