Water and Biomolecules

1   The N-methylacetamide (NMA) dimer

In this section, we look at different models for peptide chains and for the hydrogen bonds between peptide chains. Before polypeptides were computationally accessible, the NMA dimer has been a popular model for a small peptide and for the peptide hydrogen bond. As the calculations take too long time to be completed in the exercise and as you already know how to calculate things, you are just given the output files.

  • What do amino acids look like, and how is a peptide chain composed of different amino acids? What does a β-sheet look like? If you do not remember this, you can refresh this with the supplementary material.

  • A relaxation of a pair of NMA dimers has been made. The results are in:

    NMA-dimer.nc NMA-dimer.txt NMA-dimer.traj

    A relaxation of a single NMA molecule has been made. The results are in:

    NMA.nc NMA.txt NMA.traj

Copy these files to your area. Look at these calculations with plottrajectory and understand the structure of the NMA molecule (see also the article, which has been handed out).
  • In order to measure bond lengths and angles in a simple way (instead of operating directly with the positions from dacapo), one can save the configuration in a .pdb file. This is done with the following python commands:

    >>> from Dacapo import Dacapo
    >>> atoms = Dacapo.ReadAtoms('outfile.nc')
    >>> from ASE.IO.PDB import WritePDB
    >>> WritePDB('plot.pdb', atoms)
    

    and the configuration is saved in the file plot.pdb. This file can be opened with Rasmol by typing rasmol plot.pdb at the command line. A Rasmol window pops up and a Rasmol command line appears. By clicking on the atoms, their number and type are printed on the command line. By typing pick distance in the command line and clicking onto two atoms, their distance appears. One can measure the angles by typing pick angle and then clicking on the three atoms, which span the angle of interest.

  • Do you think that the NMA dimer is a good model for a peptide chain or a β-sheet? Which features of a peptide chain does it include and which are not modeled well?

  • Calculate the interaction energy between the two NMA molecules. Which hydrogen bond is mainly responsible for the attractive interactions? Would you expect the interactions of a real β-sheet to be stronger or weaker than the interaction between two NMA molecules? Calculate the interaction energy for PW91 (the functional used in the calculation) and also for PBE and RPBE (just take the non-selfconsistent energies). Try to comment on the differences.

  • Two static calculations have been carried out, where each of the two NMA molecules forming the dimer are frozen: The structures are not relaxed, but they are kept in the position, in which the molecule forms the dimer. The results are contained in the files:

    NMA_freeze1.nc NMA_freeze1.txt

    for the first NMA molecule and in:

    NMA_freeze2.nc NMA_freeze2.txt

    for the second NMA molecule. Copy these files to your area.

  • Now we are interested in, how the electronic density is changed, when two NMA molecules interact, i.e. how the two molecules are polarized. For this we plot the difference between the electron density of the NMA dimer and the sum of the densities of the two single NMA molecules alone. This is done by the script:

    diff_density_plot.py

    Try to understand, what this script does. Then run it interactively in a python prompt.

  • In the plot, blue and yellow areas show up. Can you understand from the script, which color stands for electron density depletion and which for electron density addition?

  • From the plot one can see which hydrogen bond is mainly responsible for the interaction. Which one is it?

2   Parallel and antiparallel two-stranded beta-sheets

A realistic extended model of β-sheets is only marginally accessible to DFT or other quantum-mechanical calculations. Therefore, calculations are often done with force fields, which are optimized to reproduce experimental β-sheet geometries. Two β-sheet structures, parallel and antiparallel have been calculated and are in the files:

gly-gly10_par.pdb gly-gly10_anti.pdb

for chains consisting of glycine and:

ala-ala10_par.pdb ala-ala10_anti.pdb

for chains consisting of alanine. Look at the structures and try to understand them. How have the peptide chains been terminated?

  • The most important hydrogen bond in β-sheets is the N-H...O=C hydrogen bond. How well do these bonds fit together in the parallel and antiparallel β-sheets, where do they fit better? Which type of β-sheet does the NMA dimer resemble?
  • Measure the bond lengths of the bonds in the peptide chain (C-N and C-C bonds) and the C=O bond. From other organic compounds the following approximate bond lengths are known
    • C=O double bond: 1.20 Å
    • C-C single bond: 1.54 Å
    • C-N single bond: 1.47 Å
    • C=N double bond: 1.27 Å

Which bonds in the peptide chain have single bond character and which have double bond character? What does that mean for the flexibility of the peptide chain?

  • Measure the bond geometries (lengths and angles) of the N-H...O=C hydrogen bond both in the NMA dimer and in the alanine β-sheets. How well do they agree?
  • What effect does water as a solvent have on the stability of β-sheets. Which type of amino acids would you expect to find in a β-sheet?

3   The water dimer

The file 2H2O.nc is the output from a Dacapo calculation of a water dimer. Study, and understand the script 2H2O.py that was used for the calculation (don't repeat the calculation).

  • The dimer structure has already been relaxed. What is the length of the hydrogen bond?

  • Make a density-difference plot like you did for the NMA molecule. Now you will have to do the two frozen water molecule calculations yourself. The first one can be done like this:

    atoms = Dacapo.ReadAtoms('2H2O.nc') # Restart from the old calculation
    del atoms[0:3]                      # Delete the first three atoms
    calc = atoms.GetCalculator()        # Get a hold on the calculator
    calc.SetNumberOfBands(6)            # Reduce the number of bands
    calc.SetNetCDFFile('H2Oa.nc')       # Write output to 'H2Oa.nc'
    print atoms.GetPotentialEnergy()    # Go!
  • What is the binding energy of the dimer? Compare with the binding energy of the NMA dimer.