Eplosions in proteins...

Main Table of Contents

Wed 27 Feb 2002

Explosions in proteins...

High energy X-rays damage matter, mainly through the photoelectric effect. A photon is absorbed by an atom and a photo electron is emitted from the inner (K) shell of the atom. The remaining hollow ion will relax through Auger decay, in which an electron from a higher level falls down and a further electron is emitted with an energy around 250 eV for atoms in the first and second rows of the periodic table. The basic atomic physics for these processes is well understood and cross-sections for photon scattering have been tabulated. Furthermore, Auger line-widths have been measured that give K-hole life times, e.g. for carbon of 11 fs. This means that it takes on average 11 fs before the Auger decay. For lighter atoms it takes longer than for heavier atoms, (mainly) due to smaller charge on the nuclues.

During the photoionization process the photoelectron may interact with electrons in the valence shells. For elements of biological significance (C, N, O, S) this may lead to so called shake-on and shake-off effects in which a further electron of low energy (10-100 eV) is emitted. Semi-empirical quantum calculations were used to estimate that such shake-off ionization occurs in 10 - 30% of the photoelectric events. These processes will be neglected here.

The advent of X-ray free electron lasers (FEL) will within a few years enable a whole range of new experiments in physics, chemistry and biology. In the latter category we envision performing structural studies on large biomolecules, biomolecular aggregates or nanocrystals. Molecular dynamics simulations of a protein molecule in a FEL beam are encouraging as to the feasibility of such experiments. (Neutze et al., 2000). See the figure below for an indication of what happens to a lysozyme molecule in the beam. The molecule explodes, but only after the 2 fs X-ray pulse is over. This means that there is some hope at least that we can recover structural information from the sample (during those same 2 fs) before the molecule explodes.


The molecular dynamics algorithm was modified for the explosion simulations. At each time step we perform the following additional steps.

  1. There is a certain probability PPhoto(t) that a photon interacts with an atom. PPhoto is dependent on time (pulse shape), pulse intensity, photon energy and atom type (different atoms have different cross sections). If we know PPhoto(t) we draw a random number R, and if R < PPhoto(t) we model an atom-photon interaction:
  2. There is a certain probability PAuger(t) that a K-hole will experience an Auger decay. PAuger(t) depends on the atom type. For all the atoms with a K-hole we again draw a random number R, and if R < PAuger(t) we model an Auger event:

Operation Terminate Protein

In the following assignment you will need to work without strict guidelines. Note that all GROMACS programs have a flag -h that gives helpful information about what the program does, and what the options are. Start by making a special directory for this assignment. Then download a protein structure from the protein databank. Try to find a protein with 300-500 residues, without cofactors and with a globular shape. Let's assume the protein structre is downloaded in a file called 1abc.pdb (note that the .pdb extension to the file name is obligatory). Now we have to generate a GROMACS topology using the pdb2gmx program (pdb2gmx -f 1abc). Select 0 as choice of force field. If you encounter problems at this stage due to cofactors (unknown molecules/residues) then fetch a different protein. If all is well you now have a file conf.gro containing the coordinates, and a file topol.top containing the molecular topology.

The simulations we want to do here are an extension to what was treated by Neutze et al., 2000. More specifically we are going to embed the protein in water layers of different thickness, and see how that affects the explosion. First we put the protein in the center of a big box: editconf -f conf.gro -o out.pdb -box 10 10 10 -center 5 5 5. Embedding in a shell of water is done by the genbox program: genbox -cp out -cs spc216 -shell XXX -o confXXX.pdb, where XXX is the thickness of the water layer. For these simulations try XXX = 0, 0.3, 0.6, 0.9, 1.2 and 1.5 nm (Hint 1: make subdirectories for each of these models. Hint 2: you don't need to run genbox for a water layer of 0 molecules.). You can view the resulting confXXX.pdb using rasmol. Verify that the molecule is covered with a water layer (using rasmol) and that no strange periodic boundary condition effects have occurred. Then copy the topol.top to each subdirectory and adjust the total number of water (SOL) molecules (genbox tells you what the number of solvent molecules is, you may have to add a line SOL N at the end of the topology file, where N is the number of solvent molecules).

Now copy the file em.mdp from the speptide subdirectory and modify it to have the line:

pbc = no
Copy the modified file to each of the explosion directories, and create an input file for an energy minimization (grompp -v -f em -c confXXX -o em). Then run the minimization: mdrun -v -s em -c after_em.pdb. Check the resulting structure file with rasmol again. If the energy has converged sufficiently we can now use these structures for explosion simulations.

A Perl script has been created to facilitate running the explosion simulations for different pulse widths and intensities. This script you can copy from this link to your working directory. If you have the correct files (after_em.pdb, topol.top) you can just execute the script: nohup ./explode.perl >& explode.log &. Now a series of 18 simulations with six different pulse widths (2, 5, 10, 20, 50, 100 fs) and three intensities (1012, 3x1012 and 1013) will be performed, and the results are located in subdirectories in a subdirectory called SIMS.

From these simulations you will first analyze the results as printed directly in e.g. ionize.xvg. You should see something like this:

For some of the simulations (in particular the longer ones) it is instructive to watch the trajectory (traj.trr) using ngmx. Then for all simulations do the following analysis:

  1. Root mean square deviation of C-alpha with respect to the starting structure
  2. Radius of gyration of the protein
  3. Radius of gyration of the whole system (protein + water)
  4. Use g_energy to extract the Potential energy, the Kinetic energy and the Total energy
For all the analyses we want to know what is the final number, i.e. at the end of the simulation. Then make a plot of these properties as a function of water shell thickness. If you e.g. consider that we get the following results for the radius of gyration of the protein:
0.0   1.5
0.3   1.5
0.6   1.4
0.9   1.35
1.2   1.33
1.5   1.32
then you make a text file containing that information, and load it in to xmgrace to view it (you can also print it). Your report should contain an analysis of all these properties as a function of water shell radius, along with your other observations and conclusions.

To make your analysis easier, a small Perl script has been created. You can copy it (by shift-clicking) here. Don't forget to make it executable after downloading (chmod +x analyse.perl). Go to the directory above all the imax* directories, and execute it. It will do all the four analyses above in each directory. Please embrace and extend the script.

As a final assignment (if there is time), you can select one of the trajectories, convert it to xtc using trjconv -f traj.trr -o traj.xtc and then visualize it using the program dino. This program can make glossy pictures, and what's more, it can make mpeg animations. Try reading the dino documentation in order to find out how to convert your simulation of choice into an mpeg animation. The nicest ones will be used on our website...