The bond energy term corresponds to the stretching and compressing of the length of a bond. In most energy functions this term reduces bonds to simple harmonic oscillators, yielding a quadratic equation: E bonds = K b ⁢b− b 0 2 E bonds K b b b 0 2 where KbKb is an empirically determined constant that depends on the atom types, bb is the current bond's length, and b0b0 is the bonds length in equilibrium, which again depends on the atom types. In this case you can think of the bond as a spring, it has an equilibrium length that it wants to remain at. If the bond length varies from the equilibrium length, the energy increases.

The bond angle energy corresponds to changes in the angle between bonds. As with bond length, the bond angles have an equilibrium value, and any deviation increases the potential energy. Once again this can be modeled by a simple quadratic term. E angles = K θ ⁢θ− θ 0 2 E angles K θ θ θ 0 2 where KθKθ is an empirically determined constant, θθ is the current bond angle, and θ0θ0 is the equilibrium angle.

Torsions are created by series of three bonds, and consist of rotations of the bonds on either end with respect to the axis of the middle bond. In molecular structure certain torsional angles are preferred over others and the energy function reflects this. Usually it is described by a Fourier series expansion. The simplest being a single term:

Figure 4: A typical torsional energy term for a potential function.

A more complicated three term expansion can also be used:

Figure 5: A more complicated torsional energy term for a potential function.

Strictly speaking, Van der Waals interactions are weak attractive interactions between atoms at an ideal separation from each other. The atoms transiently induce each other's electron distribution into complementary dipoles, allowing a weak electrostatic attraction between them. In molecular potential fields, Van der Waals attractions are usually combined with steric clash (extremely high energies due to overlapping atoms) in a Lennard-Jones potential, such as this Lennard-Jones 12-6 function:

Figure 6: A typical Lennard-Jones 12-6 potential.

Electrostatic interactions are usually computed using some variant of Coulomb's Law, which assumes that atoms behave as point charges located at their centers. A typical Coulombic term looks like this:

Figure 7: Electrostatic energy is computed using a version of Coulomb's Law.

The dielectric constant is a function of the medium through which the two charges interact. The difference between the dielectric constant of water and that of pure protein is substantial, so some models attempt to take it into account. One of the simplest assumes that the farther apart two charges are, the more likely they are to have water between them. This is called a distance-dependent dielectric, because it scales with the distance between the atoms involved.

While all of the previous terms are almost always included in energy functions, there are a handful of other terms that are common, but not present in every function. These include hydrogen bonding, solvation and cross terms.

Hydrogen bonds (which are not true bonds in the strict, electron-sharing sense) are unusually strong electrostatic interactions, usually between a hydrogen atom and an electronegative atom such as oxygen or nitrogen. They play an important role in determining and maintaining the structure of biomolecules including proteins and nucleic acids. Some energy functions account for hydrogen bonding in the electrostatic term. Other functions include a separate hydrogen bonding term which is most often a Lennard-Jones-like 12-10 potential:

Figure 8: A hydrogen-bonding 12-10 term for potential functions.

The solvent that a molecule is in can have a large effect on how it moves. Explicitly representing solvent molecules, however, is a computational cost that most methods try to avoid. Usually the solvent model is separate from the energy function. There are several different ways of approximating solvent interactions including the Generalized Born Model and the Poisson-Boltzmann method. Most force fields do not have an explicit solvent term.

Other terms that describe the interaction between bonds and angles, angles and torsions and so on are included in some force fields. For example to model the interaction between bonds and angles:

Figure 9: A potential energy term depending on both bond lengths and angles.

The research group of Nancy Amato has been working on roadmap-based methods to study the process of protein folding [1][2][3][4]. They start with the known native structure of a protein, and incrementally find conformations more and more different from the native state, and build a roadmap using these conformations. The goal is to find a large ensemble of pathways between the native structure and unfolded structures, and to study these pathways and their properties. In their work, the degrees of freedom of the protein are assumed to be the φ and ψ backbone dihedral angles of each amino acid residue. The side chains are assumed to be rigidly attached to the backbone. In their initial work, they generated new conformations for their roadmaps by adjusting the backbone dihedral angles in the folded conformation randomly with various standard deviations. This approach worked well for very small proteins, but did not scale well.

A later sampling approach was based on counting native contacts. For the purposes of their method, a native contact was defined as any two alpha-carbons within 7 Å of each other in the folded state of the protein. A new conformation is generated at each step of the sampling phase of the roadmap construction by randomly perturbing some existing sample. The resulting structure was accepted with a probability as follows:

Figure 11: The acceptance criterion for newly sampled conformations.

Once the samples are generated, an attempt is made to connect the k nearest neighbors of each node to the node itself. A series of conformations on the line connecting the two nodes are tested, and if their energy is below a threshold, the edge is included in the roadmap. The weight of the edge depends on the sequence of energies of the conformations computed along the connecting line. The probability of moving from the (i-1)th structure to the ith along the line is given by:

Figure 12

Once the roadmap is computed, the shortest paths between structures can be found using Djikstra's algorithm, and the folding paths can be extracted and studied. In particular, the order of secondary structure formation can be predicted by a consensus method. The order of secondary structure formation is determined for all paths in the roadmap from unfolded to folded structures, and the most common order is predicted as the true order of formation, which is a coarse, high-level expression of the folding mechanism. On a set of proteins used to test their method, the predicted formation order matched laboratory-determined formation order in all cases where it was available. Because of the coarseness of this notion of the folding mechanism, a statistical analysis of all pathways makes sense.

In their most recent work [4], these researchers have refined the method by which new structures are generated in the sampling phase of roadmap construction. This method is based on rigidity analysis. For each structure, each degree of freedom is determined to be independently flexible, dependently flexible, or rigid, using an algorithm called the Pebble Game [13]. Independently flexible degrees of freedom rotate with no deterministic effect on others. Dependently flexible degrees of freedom force other degrees of freedom to change in a set way. Rigid degrees of freedom are essentially locked in place unless the constraints change. In perturbing an existing structure to generate a new sample, degrees of freedom are perturbed with a strong bias towards perturbing independently flexible degrees of freedom and against perturbing rigid degrees of freedom. Because the new structures are generated by physically realistic motions, it is expected that they will generally be lower energy and more representative of real structures than if they were generated by completely random perturbation of the degrees of freedom.

In practice, rigidity sampling appears to allow construction of high-quality roadmaps with many fewer samples than were necessary without it. It thus improves the overall efficiency of calculating protein behavior using this roadmap method.

Numerous variants of MD and MC have been developed in an effort to speed up the process or focus the simulations on particular motions of interest. One method, called the Stochastic Roadmap Simulation (SRS) [5] [6] [7] [8] uses a PRM-like structure to approximate a large number of simultaneous MC simulations very rapidly, allowing the analysis of an ensemble of trajectories. This method followed very directly from the first roadmap studies of molecular properties by Singh and Latombe.

The SRS method proceeds as follows:

The transition probabilities are defined as they are to be consistent with a Boltzmann-like distribution of energies, and therefore with standard Monte Carlo simulation probabilities. The authors demonstrated that each continuous path in the roadmap may be interpreted as a Monte Carlo simulation, and that, if a very large number of samples and edges are made, the aggregate behavior of these Monte Carlo simulations can be analyzed to estimate properties of the protein such as folding rates and transition states. Essentially, SRS is a way to generate large amounts of Monte Carlo simulation data in a short time. The developers or this method have provided a proof that, for a sufficiently large SRS and a sufficiently long Monte Carlo simulation, the distribution of conformations is expected to be equal.

To study protein folding using SRS, we observe that some set of nodes in the roadmap represent conformations in or very close to the folded state (native structure). We will refer to this set of nodes as F. For every node in the roadmap, we can compute an expected number of state transitions (or Monte Carlo steps) to go from that state to a node in F, with the base case that any node in F is defined to be at distance 0 from F. Given a precomputed SRS, we can compute this statistic for each node as follows:

Figure 15: The expected number of transitions to reach a node in the folded state starting from node i.

We can also define a set of nodes representing conformations close to the stable denatured state of the protein as the unfolded state, U. Given both of the sets U and F, we can define a quantity called the transmission coefficient, τ, for each node. The transmission coefficient expresses the probability that a structure at a particular node will proceed to a state in F before it reaches a state in U--in other words, it is the probability that a given structure will fold before it unfolds. This is often called the folding probability, or Pfold, in more recent research. The quantity, τ, is calculated for each node using the following relation:

Figure 16: The folding probability for a node i. This is the probability than a simulation starting at node i reaches a folded state before reaching an unfolded state.

Markovian State Models (MSM) [9] [10] are roadmaps constructed by running many molecular simulations (Monte Carlo and molecular dynamics) and merging the trajectories. The method starts with a simulation that runs from the folded to unfolded state. It then picks a structure at random (call it c) from this trajectory and run a new simulation. If this new simulation reaches the unfolded state, then the next trajectory from which we will pick a structure will consist of the old trajectory from the folded state to c, and the new trajectory from c to the unfolded state. If the new trajectory reaches the folded state, we do the opposite. If neither happens in a reasonable time, we reject the new trajectory and start over. This is repeated a set number of times, and each time a trajectory is accepted, all states from the new part of the trajectory are added to the growing roadmap as nodes, and each transition from the trajectory is added as an edge. When it is added, each edge is labeled with a transition probability of 1 and a transition time equal to the timestep of the simulation.

The goal of this method is to move roadmap methods closer to MD and MC simulations, and in particular to incorporate a notion of time, which follows from the use of simulation techniques in the sampling of new conformations.

Once all of the simulations have been run, nodes that are within some cutoff distance of each other by some similarity metric must be merged. To merge two nodes, one node is removed, and its edges are transferred to the other node. If this results in two edges between the same pair of nodes, the transition probabilities and times are defined as follows:

Figure 17: Expressions for new transition probabilities and transition times when merging nodes in constructing a MSM.

An important fact of all roadmap methods that attempt to extrapolate properties of the entire protein folding landscape is that there is inherent sampling error. The energy landscape of a protein is a continuous function, which roadmap methods attempt to approximate through discrete sampling. The researchers who developed the MSM method also developed a method to estimate the error of the folding rates estimated based on MSMs [10]. While a complete description is beyond the scope of this module, the details are available in the 2005 paper by Singhal and Pande linked in the Recommended Reading section below. The error analysis allows them to generate a probability distribution for the folding times for each node in the MSM. Useful in its own right because it gives us an idea of how confident we can be in folding times generated by a given MSM, this analysis is especially useful for focusing sampling during the generation of an MSM.

The variance of the distribution of the folding time for each node provides an estimate of the error. If at each stage of simulation instead of choosing a node at random to start the next simulation, we select the node with the greatest contribution to our estimate of the error in folding time, we effectively focus our efforts where they will decrease the error most. In this way, MSMs with less overall error may be generated with using simulations.