Refinement and stability of protein structures: an application of MD

The optimal situation would be if a computer simulation of the folding process could be able to predict the structural properties of proteins. Unfortunately, an exclusively molecular dynamics (MD) based method works only for some small peptides having stable structures. The origins of this restricted applicability are two-folded. On the one hand the available computers and methods are unable to follow the folding process for an appropriately long time at a proper level of accuracy. On the other hand the details of the structure formation are not completely known experimentally: e.g. the details of interactions with other molecules, like chaperons.

Therefore we often use initial structures obtained either from experiments (like X-ray or NMR) or from theoretical predictions (like homology modelling). However, these structures frequently need refinements, because the experimental conditions are far from the physiological ones or the quality of the homology model is not good enough for the further investigations etc.

In this example MD simulation of a mini protein (Trp-cage) is used for refining the structure as well as to investigate its stability. Root-mean-squared-deviation (RMSD) and Root-mean-squared-fluctuations (RMSF) of atomic coordinates are used as mathematical tools in our study.

3.1. Comparing to a reference structure

Let us suppose that we use an experimental or theoretically derived structure of a protein as a starting point in an MD simulation (with a proper solvent, temperature and appropriate co-soluted ions). During the simulation this structure goes through smaller or larger rearrangements. If the rearrangement is small, then the initial structure is close to the equilibrium structure of the protein at the simulated conditions (supposing that the force field parameterization and other parameters used provide good approximations of the physiological conditions). The large rearrangement shows the opposite. For a quantitative study we need a mathematical measure of the difference between the two structures.

3.2. RMSD,least square fitting

Let us denote the coordinates of the i-th atom of our protein r _i(t₀) and r _i(t) at the beginning (t₀) and at a time point t of the simulation. The differences of these two structures can be characterized with their root-mean-square deviations defined by the following formula

(13.1)

The structure of proteins is dominantly determined by the heavy (non-hydrogen atoms). The positions of H-s are less important from this point of view than other atoms in the characterization of structural similarity. Mathematically this can be formulated using the mass weighted RMSD (MW-RMSD) defined as

(13.2)

If we calculate the RMSD for each residue separately we can follow the distance of a residue from its reference position along the MD trajectory and identify the time and places of the largest structural changes, too:

(13.3)

(13.4)

(13.5)

The simplest sign of the stability of a structure is that it does not change significantly during the MD simulation (equilibrium structure). But even in this case the geometry fluctuates around an average value. The extent of these fluctuations can be characterized by its average squared deviation of atomic positions from their average values, which is also called variance or mean square fluctuation (σ, MSF):

(13.6)

(13.7)

3.4. MD investigation of Trp-cage miniprotein

To demonstrate the ideas described above, MD simulations were performed on the Tryptophane-cage (Trp-cage) peptide (it is also called TC5B), using explicit water molecules to model the surroundings of the protein. Trp-cage is a 20-residue-long peptide with the sequence NLYIQ WLKDG GPSSG RPPPS, which is often called miniprotein, since it is one of the smallest polypeptides that possesses a well defined secondary structure [1]. The three-dimensional structural features of this peptide are well characterized by several experimental methods and molecular modelling techniques. The most important structure stabilizing factors of Trp-cage are the hydrophobic stacking of the aromatic rings of Tyr3 and Trp6 amino acids, and the salt bridge formed between the Asp9 and Arg16 residues. Its small size, temperature-sensitive structure and the simultaneous appearance of two important stabilizing interactions make this system an ideal protein model for MD studies.

The Amber _99SB-ILDN [2] force field and the TIP3P [3] water model were used in our NPT ensemble MD simulations using the GROMACS [4] molecular dynamics package.

Figure 13.5. Structure of Trp-cage (tc5b) mini-protein. The secondary structure is represented using the “new cartoon” style of VMD [5] program

Figure 13.6. A: values of Trp-cage miniprotein. The average value of 1.89 Å for the last 46.6 ns is shown by the red line. B: The experimental (red) and average (yellow) structures of Trp-cage C: Structural fluctuations around the average structure: snapshots were taken from the trajectory

Figure 13.7. d^R_{aMW - RMSD} values of Trp-cage miniprotein for the backbone atoms and for all atoms of each residues, separately. The experimental structure was used as reference.

Figure 13.8. RMSF values of Trp-cage miniprotein for each residues

Figure 13.9. A short animation of the structural fluctuations around the average (yellow) structure.

4. Binding affinity estimation

The free-energy change associated with the binding of a ligand to a protein can be estimated with the Linear Interaction Energy (LIE) method [] (see Chapter 9). This method appears to be a good compromise between accuracy and computational efficiency but it does not perform equally well for all systems. The first step in its application is the determination of some system dependent parameters by using experimental binding free energy data of a set of related compounds all bound to the same protein. This fitting does not only provides us with the parameters required but is also gives us information on how well LIE works for our system.

The application presented below estimates the binding free energy between prolyl oligopeptidase (POP) and some of its ligands. It will also show that the calculations give insight into the details of the binding and the relative importance of electrostatic versus van der Waals interactions. These pieces of information may be exploited in ligand design.

POP is a serine protease that cleaves the peptide bond on the carboxy side of proline residues. A series of its inhibitors is shown in

Figure 13.10. Structure if prolyl oligopeptidase inhibitors studied. Reprinted with permission from J. Med. Chem., 51, 7514–7522, (2008). Copyright 2008 American Chemical Society.

Figure 13.11. Binding of the P1-P2 moiety of the inhibitors in the crystal structures. The ribbon model of the protein is colored gray, while the inhibitor molecules and POP binding sites of the P1-P2 moieties are magenta, red and green for the POP-6, POP-8 and POP-11 complexes, respectively. Hydrogen bonds are shown as shaded lines. Reprinted with permission from J. Med. Chem., 51, 7514–7522, (2008). Copyright 2008 American Chemical Society.

(13.8)