Computational Chemistry Courses

Normal Mode Analysis (Tama and Sanejouand, Protein Eng. 2001; Dobbins et al, PNAS 2008)

Linh — Fri, 17 Apr 2009 18:33:00 GMT

Proteins with many transitional conformations when changing from one conformation to another usually involve relative motions of domains. A good understanding of these movements will help the study of protein functions such as catalysis and regulations. In these papers, a theoretical method for studying collective motions in proteins, normal mode analysis (NMA) is used. With this tool, conformational changes can be expressed in term of superposition of collective variables, ie normal mode coordinates.

Normal Mode Theory:

The displacements of atomic coordinate i, r_i(t) with near by stationary point of the potential energy surface can be calculated using Goldstein equation (1950). Within this approach, a simplified potential energy function proposed by Tirion (1996) is used. This simplified potential energy function is designed such that any configuration of any system is a minimum of the function. With this method, NMA no longer require energy minimization step and thus, cut off CPU time. Another two techniques of the NMA in this study are all the normal mode calculation used a cutoff of 8 Å which has been used by Bahar et al (1997) and only C_α atoms were taken into account. This model allowed the study of backbone motions, which characterizing low frequency normal modes of large proteins, to be achieved with small amount of CPU time and lower computational expense. The simple potentials and models have been used in previous studies, and the results strongly suggested that simple potentials yielded low frequency normal mode as accurate as those with standard NMA.

In order to compare with the experimental method, the authors looked at two qualities, the overlap I_j and the correlation coefficient C_j where I_j is the overlap between delta_vector_r = {Δr₁, ..., Δr_i, ..., Δr_3N}, the conformation change in crystal structures and the jth normal mode of the protein; C_j is the similarity of the patterns of atomic displacements in the conformational change and in mode j. Equations for quantitatively calculations of C_j and I_j are given in the paper (p2-3). Note that delta_vector_r is calculated for pairs of crystallographic conformations of “open” and “closed” form after both conformations were superposed. The overlap and correlation are considered as functions of degree of collectivity (к) of the conformational change where degree of collectivity is a measurement of how collective a motion of the protein is (Bruschweiler, 1995). Overlap of the most involved modes in conformational changes calculated using simplified model are found to be almost equivalent as those obtained using standard NMA and to have smaller RMSD 1.2 – 1.9 Å comparing to 2-3 Å obtained in standard approach. The results also suggested that NMA performed better in the open form than closed (refer to table III), better with highly collective motions (refer to table IV) and better with more localized motions (refer to the case of tryglyceride lipase in table IV and V). In conclusion, a single normal mode can carry a lot of information about conformational transitions of proteins. The single mode with best correlation and best overlap is more likely not the lowest frequency mode but often one of the three lowest frequency modes.

Carrying on the NMA to study conformational change of protein, the work in the second paper performed NMA calculations to investigate the extend to which thermal motions of proteins could provide the general nature of conformational change upon protein-protein docking.

A set of 20 proteins with flexibility which have been observed to have large conformational change under forming complexation by Tama et al, is used in this study. Their result showed high agreement with the experiment in term of average collectivity, 0.40 for observed and 0.38 for predicted. They also stated that the RMSD of binding sites are found to be 1.2 times larger than those over the entire protein, which suggested that upon complexation, the protein conformational change has highly flexible regions around binding site. Further investigation of mobilities over the binding site for both observed and predicted using NMA were performed. Splitting the binding site into the core and the periphery residues, they found that the observed motion of the peripheral residues is more than twice large than that of the core residues. This confirm the hypothesis that binding site has regions of highly flexibility along with rigid regions upon docking (refer to figure3). There were no obvious trend for the overlap. They were observed to be different for different frequency modes. However, their results are in agreement with previous study that large conformational change of proteins usually have low frequency modes but not the lowest frequency mode. The protein size are also examined in this paper. They proposed that the larger the protein size the more collective motions it has and thus might also be used to predict conformational change. However the size is less reliable than the normal mode frequency.

Overall, this particular study showed that almost all proteins undergo large conformational change has thermal motions in isolation and required some intrinsic flexibility. In their prediction using NMA, only one third of the proteins have the direction and location of motion among the lowest modes is similar to the observed conformational transition. They suggested this because of the bound condition which might altered structural changes. In my opinion, this is one of the limitations of this work, along with unclear trends for overlap. In conclusion, the significant contribution of this study is that it suggested a model to study the conformational change during molecular recognition will be desirable in future work.

References:

Tama and Sanejouand, Protein Eng., 2001

Dobbins, Lesk, and Sternberg, PNAS, 2008

λ-dynamics method and applications (Kong and Brooks, JCP 1996; Khandogin and Brooks, Biophys J 2005)

StevenWang — Thu, 16 Apr 2009 14:52:00 GMT

λ-dynamics is a new method developed on the basis of the widely used free energy perturbation
method (FEP) and the umbrella sampling method. It differs from the FEP method by introducing
multiple parameters λi, instead of only one λ, which are assigned to different types of energy
components. In order to surmount the possible energy barriers and thus get more samples,
special biasing potentials are used based on the idea of umbrella sampling. The forms of these
biasing potentials are quite innovative, which look like kinetic energy and potential energy for
factious particles where λi, instead of r, serve as the coordinates and a set of fake masses of
particles are introduced. The benefit from introducing more λ is that more reaction pathways are
accessible during the simulation. On the other hand, the umbrella potential term is useful in
flattening the energy curve of the reaction pathway, preventing the boundary crossing of λi (e.g.
λi < 0 or λi > 1) and limiting the range of {λi}. Additionally, the fake masses make the control of
λi more flexible. In a word, these new variables and terms make the simulation more adjustable
and easier to control for different purposes. The results of this new method are at least as
accurate as those by the single λ FEP method, which are verified by simulations of relative free
energies of hydration and relative binding affinity simulations.

Constant pH molecular dynamics (CPHMD) is one of the applications of the λ-dynamics
approach. In CPHMD method, n λi are assigned to n titrating residues, which varies from 0 to 1
when the titrating residue change from protonated form to deprotonated form and it takes the
mathematical form as an sine function. A two-dimensional λ-dynamics method is used here for
the purpose of take the proton tautomerism effect into consideration. The two dimensions are
due to the two reaction coordinates: one is for the deprotonation process (λ); the other is for the
tautomeric effect (x). These two coordinates are used throughout the potential energy
calculation via λ-dynamics method. When calculating van der Waals energy, pairwise
interactions are divided into two classes: interactions between titrating residues and other
resides and those between two titrating residues. For coulomb and GB electrostatic energies,
atomic partial charges on titrating residues are changed along the reaction coordinates. Biasing
potentials are also used to facilitate the simulation. One is like harmonic umbrella potential and
the other one originates from the calculation of free energy based on experimental pKa value.
In CPHMD, the unprotonated form has λ > 0.9, while that of the protonated form is smaller than
0.1. In both cases, the x value is either larger than 0 or smaller than 1 by a maximum magnitude
of 0.1, in other words, the tautomer state must be sufficiently pure. The enhance the sampling
near the two ends of the reaction coordinates, a barrier potential is introduced. This potential
function is symmetric considering the middle point of the reaction coordinate. The result is that
the energy is lower when the titrating residues are in pure states than that of the titrating
residues in mixed states and therefore more pure sate residues are sampled.

In my opinion, the CPHMD method makes a best use of the advantage of the λ-dynamics
method, making the titration simulation more flexible and closer to the real experiments. Thus , I
expect more repeatable and accurate data will come from this method. But there are still some
minor problems. The authors assume the van der Waals energy changes linearly or
quadratically during the change of reaction coordinates, which may not be true compared with
that in the real world. Besides, the rate of the changes of van der Waals energy is coherent with
that of electrostatic energy in this paper. It is possible that assigning different reaction
coordinates to them may produce better results, just like what is done in the first paper.

References:
Kong and Brooks III.1996. λ−dynamics: A new approach to free energy
calculations. J.Chem. Phys. 105:2414-2423

Khandogin and Brooks III. 2005. Constant pH Molecular Dynamics with Proton
Tautermerism. Biophysical Journal. 89:141-157

Challenges in Theoretical Chemistry - Intermolecular Potentials (Stone, Science, 2008)

StevenWang — Thu, 16 Apr 2009 02:49:00 GMT

Main points: This paper concerns about the basics, as well as the capabilities and limitations of the intermolecular potential calculation. As a novice, I like the basic part better. Because the small molecules cases are easy to cope with, the author cares more about the intermolecular potential calculation between two large molecules. In this situation, the author claims, intermolecular potential calculation is transformed into pairwise atom-atom potential calculation, which is the sum of all the potentials between any of the atoms in molecule A and another one in molecule B. This can cause a double count problem but the author does not mention too much about it. The calculation of the potential is through a general-form function: U {sub a & b}, whose variables are interatomic distance R and relative orientation Omega between the two atoms involved. For the Omega (relative orientation), I think this the result of concerning the electron distribution in the electron cloud of each atoms or the partial charge distribution, which is discussed later in the paper.

There are a total of five terms contributing to the potential function: 1) the electrostatic term; 2) the exchange-repulsion term or van de Waals repulsion; 3) the dispersion term or van de Waals attraction; 4) induction term; 5) charge transfer term. (1) For the electrostatic term, the author reviewed several methods that involve charges, dipoles, multipoles, etc, instead of using point charges. And, truly, these methods can interpret the electrostatic term better, especially in the hydrogen-bonded systems. (2) For the exchange-repulsion term, according to the author, the Born-Mayer form works better than the R {superscript -12} term in the Lennard-Jones potential. (3) As for the dispersion term or van de Waals attraction, it is calculated through a series with a introduction of a correction term (a damping function) for R near zero. (3) The induction term and charge transfer term are the most troublesome ones, from the perspective of the author. The induction term comes from the interaction between an electric-field-giving atom and the consequent induced dipole of another atom. The charge transfer term is caused by the donating of electron density of one atom to the acceptor molecular orbital (MO) of another atom (maybe I should say the MO of another molecule), which will give rise to charges in electron configuration. There are two ways of dealing with this induction term: inexplicit induction, by modifying the properties of each molecule and explicit induction, by assigning isotropic dipole-dipole polarizabilities to atoms and repeating solving induction term at each step during the optimization or simulation.

Up until now, it is already practical to accurately calculate the interatomic potentials for small molecules (using, say, ab nitio method) and molecules made up of 20-30 atoms (using, say, symmetry-adapted perturbation theory based on DFT). For macromolecules, the simplest way is employing atomic point charges, Lennard-Jones or Born-Mayer repulsions, and isotropic atom-atom dispersion, while other treatments are also available. However, there is currently no definitive all-purpose force field.

Two limitations are associated with intermolecular potential calculation. One is the difficulty in describing flexible molecules. The other is the limited available experimental data for the use of parametrizing the force field.

Generally, this is a good review of the development of intermolecular potential calculation and it is appropriate for all level reader interested in this field.

References:

Stone, Science, 321, 787-789, 2008.

Replica-Exchange MD (Sugita and Okamoto, Chem Phys Lett, 1999; Nymeyer et al, Methods Enzymol, 2004)

Kevin — Tue, 14 Apr 2009 22:22:00 GMT

The purpose of replica-exchange molecular dynamics is to provide adequate sampling of a wide conformational space MD simulation. In MD simulations, especially for biomolecules, the energy landscape often contains large energy barriers that will trap the system in saddle points, or local minima. In order to find the absolute minimum of the system we require a method that will allow us to cross over the large energy barriers in search of deeper minima. This is where replica-exchange helps us. In replica-exchange we create several copies of a system that are identical except for their temperature. The higher temperature replicas will be better able to sample over the entire conformational space and not be trapped as easily in local minima. Each replica is searching for energy minima. After a set number of steps in the simulation the energy of each replica is compared with its nearest neighbor(s). Exchanges of temperature are then determined according to the following probability: min{1,exp(ΔH*Δβ). Thus, if a higher temperature replica is at a lower energy than the lower temperature replica, the temperatures will be exchanged and the simulation continues. If the energy is greater there is still a chance, although small, that the temperature will be exchanged. The basic rundown is this: if a high temperature replica finds a deep energy well, the temperature will be lowered for the newly found well repeatedly until the absolute minimum is reached. If an even deeper energy well in the landscape is discovered by a higher temperature replica, the process repeats itself focusing on the new energy well. This continues until, ideally, the absolute minimum of the energy landscape is found.

References:

Sugita and Okamoto, Chem Phys Lett, 1999
Nymeyer et al, Methods Enzymol, 2004

Challeges in Theoretical Chemistry - Surface Scattering Simulations (Kroes, Science, 2008)

Kevin — Tue, 14 Apr 2009 19:10:00 GMT

The goal of surface scattering simulations is to accurately model how molecules interact with surfaces. This is a daunting challenge because molecule-surface interactions can result in many different processes including: dissociation, scattering, and various excitations. Currently most molecule-surface interactions are modeled using the Born-Oppenheimer approximation which limits the reaction to one potential energy surface (PES). However, the BO approximation does not generally hold for non-adiabatic processes. Recent research has studied whether or not there are certain reactions that can be accurately described using the BO approximation. In general the BO approximation has proven fairly accurate for light molecules (i.e. H₂) but heavier molecules tend to experience other surface-interactions that affect its scattering such as phonon excitation and electron-hole excitation. These scenarios ideally should be modeled using a diabatic PES or multiple PESs which is not the case with the Born-Oppenheimer. Many challenges remain in this field including: developing an accurate method to model interactions in the electronic ground state, excited electronic states, as well as accurately describing the effects of phonon excitation. Overall the article does not discuss surface scattering as much as I expected. Rather it focused heavily on modeling phonon and electron-hole excitations. These do present a challenge to calculating accurate scattering simulations, but I expected more discussion on how scattering interactions (vibrational, rotational, and diffractive) are calculated.

Reference: Kroes, Science, 321, 794-797, 2008

Challenges in Theoretical Chemistry - Quantum Dynamic of Chemical Reactions (Clary, Science, 2008)

Linh — Fri, 10 Apr 2009 21:25:00 GMT

Quantum Dynamic (QD) is used to study chemical reactions in the gas phase.Current computational and theoretical methods are used to calculate the rate of reactions and to predict the detail of the reactions.

Methods: For reactions with high energy barrier such as those involving bond breaking or formation, QD methods need to be used.

QD procedure: Solve Schrodinger equation to get the electronic structure. Compute Potential Energy Surface (PES) using ab initio electronic structure methods. This yield reliable reaction barriers. Use these PES in the solutions of Schrodinger equation for the nuclei to predict fundamental informations of chemical reactions.

For three atoms reactions, a computer code called ABC, which involves hyperspherical coordinates, is used. This theory can be extended to full dimensions by including functions that describe molecular rotations. Agreement between this theory and experiment is remarkable. For reactions with several maxima and minima in PES, "direct dynamic" methods are sometimes used. In these methods, electronic energies are calculated directly from an ab initio code whenever they are required. However direct dynamic are difficult to aplly in QD calculations due to large number of quantum wave functions. An alternative method is used which involves solving the time dependent Schrodinger equation. This is called the "wave packet" calculations.

For four or more atoms reactions, a method is developed that allows the construction of effective Potential Surface (PS) in reduced dimensionality through a small number of accurate electronic structure calculations. These PS are used in time independent quantum scattering calculations with hyperspherical coordinates, treating only bonds broken and formed in reactions explicitly. This method is used in reaction between H atom and propane which has result agree swell with experiment. Another approach is multiconfigurational time dependent Hartree (MCTDH) which treat quantum wave functions as a product of single particles and time dependent functions. MCTDH can be used with a method for calculating the quantum fluxes in a chemical reaction to calculate the reactive rate constant of multidimensional system.

Challenges in Theoretical Chemistry - Modeling Materilals Properties (Carter, Science, 2008)

Jason — Tue, 17 Feb 2009 15:52:00 GMT

This perspective article in Science Magazine outlines the goals and challenges in attempting to model materials properties from first principles. The article is a concise review of the specialized techniques that are most successful in modeling different types of properties. Reading this article would be a good first step to find the method of choice to model a particular property of a material. Of course, as a review, the article does not discuss how the different techniques are actually implemented or why one technique is better suited for a particular problem. You must go to the references to get a deeper level of understanding for each technique mentioned.

Reference: Carter, Science, 321, 800-803, 2008.