Biotin/Streptavidin Tutorial

The streptavidin protomer is organized as an 8-stranded beta-barrel. Pairs of the barrels bind together to form symmetric dimers, pairs of which in turn interdigitate with their dyad axes coincident to form the naturally-occurring tetramer.

In this tutorial, we massage a PDB file of the tetramer, solvate the region of interest (i.e. around one of the four complexed biotins), and run some dynamics keeping the rest frozen. The equilibration takes 2.3 hours per picosecond on moderately fast CPU (Convex 3820). Allowing the whole system to move takes 4 hours/psec (the whole system would have to be solvated, thus still longer time, for this to be useful). The frozen part provides a more realistic electrostatic environment for the part that moves.

The initial biotin/streptavidin tetramer was prepared by Richard Dixon from the momomer from Brookhaven 1stp by P.C. Weber, D.H. Ohlendorf, J.J. Mendolowski, and F.R. Salemme (1992).

Related material, not necessary to understanding the tutorial:

   Miyamoto S, Kollman PA.
     What determines the strength of noncovalent association of
	ligands to proteins in aqueous solution.
     Proceedings of the National Academy of Sciences of the 
	US of A, 1993 Sep 15, V90 N18:8402-8406.

   Miyamoto S, Kollman PA.
     Absolute and relative binding free energy calculations of 
	the interaction of biotin and its analogs with streptavidin 
	using molecular dynamics free energy perturbation approaches.
     Proteins-Structure Function and Genetics, 1993 Jul, V16 N3:226-245.

Overview

• Prepare pdb
• Leap:
  • Build biotin residue template
  • Load 'frcmod' file (extra force field params)
  • Load pdb
  • Add 'cap' of water
  • Save 'parm' files for dynamics
• Carnal:
  • Figure out residues in biotin 'cap' region
• Sander:
  • Minimize and run dynamics to equilibrate
  • Check equilibration
  • Run more dynamics
• Leap:
  • Analyze the trajectory

Prepare pdb

	Hydrogen naming conventions in the given pdb file were 'wrong' -
	a not uncommon experience. After finding it tedious to correct the 
	names, I deleted all H's (about 4000) by using 'egrep -v' to exclude 
	lines matching a pattern with H in either of the 1st 2 columns of the 
	atom name:

	    % egrep -v '^.............H' given.pdb > x
	    % egrep -v '^............H' x > start.pdb

	where the '^ starts the pattern at the beginning of the line and
	the .'s are wild-card single characters. note


Leap
  Prepare Biotin residue template
       From prep file
       From pdb file
    (if this was done in a previous session,
		> loadoff btn.lib
    to reload)

  Load 'frcmod' file for biotin:

     > fmod = loadamberparams btn.frcmod

  Load prepared pdb of streptavidin/biotin complex:

     > stbt = loadpdb start.pdb
     > edit stbt

	


  Hold down the two right buttons and push forward/back to zoom in/out.
  The middle button alone rotates, the right button translates, and
  the spacebar recenters the molecule.

  Add a 'cap' of waters around the site of the 1st biotin. This is
  done by estimating a median x, y, z coordinate by eyeballing the
  coordinates of the BTN at the beginning of start.pdb:

     > solvatecap stbt WATBOX216 { 35 12 -6 } 20

	


  The number of waters in the cap varies slightly, depending on 
  machine; about 265 is to be expected on SGI; the HP used for
  this demo got 277.
  
  Save system in leap and pdb formats for future reference:

	> saveoff stbt built.lib
	> savepdb stbt stbtcap.pdb

  Save 'parm' files for dynamics & perturbation:

	> saveamberparm stbt stbtcap.top stbtcap.crd

Carnal

  For dynamics, we only want the region of interest to move -
  the rest is there to provide a more lifelike environment.

  Use a beta Carnal feature to figure out the residues around
  the biotin molecule that the water cap is on:

	% carnal < carnal_grp.in > carnal_grp.out

  The resulting list of residues (carnal_grp.out) is in the GROUP 
  format described in Appendix B of the manual.

Sander

  Minimize a little:

	% sander -O \
		-i min.in \
		-p stbtcap.top \
		-c stbtcap.crd \
		-o min.out \
		-r min.rst

  Warm gradually, using 'belly' option to restrict motion to
  region of interest as determined by Carnal, above. (The
  group defined by Carnal is modified in md0.in to include all 
  the waters.)

	% sander -O \
		-i md0.in \
		-p stbtcap.top \
		-c min.rst \
		-o md0.out \
		-x md0.crd \
		-r md0.rst

  This equilibration immediately warms to 300K; this may be ok
  in this case since a major part of the system is frozen - if 
  this were not the case, too-rapid warming could disrupt the
  structure. Similarly, if this was a constant pressure simulation
  in a periodic box of solvent, too-rapid warming could lead to 
  excessive pressure fluctuations.

  Although no such critical problems can happen here, we still 
  need to see if the system has truly equilibrated:


				Equilibration: Temperature

	


  The temperature seems fairly equilibrated, but this does not
  necessarily mean the system is fully equilibrated, especially
  because the temperature coupling we use in AMBER forces the
  temperature to the desired value. (Note: since this setup uses
  an effectively infinite cutoff, once it is equilibrated the
  temperature scaling could be turned off; see NTT in the 
  Sander manual.)

  Looking at the overall energy:

				Equilibration: Energy

	


  Here we see a definite trend to lower energy that presumably
  has not completed. Consider the main components, kinetic and
  potential energies:

				Equilibration: Kinetic Energy Component

	


				Equilibration: Potential Energy Component

	


  Thus after the initial warming period, kinetic energy is relatively
  constant (it scales with temperature so is not really worth checking,
  in fact), but the system is gradually relaxing to a lower-energy 
  conformation. Further equilibration is called for.  Note: in a
  periodic system, we would also want to check pressure fluctuations and
  density (see methane. 

  We restart from the saved coordinates of the previous run, using the 
  same conditions but setting the appropriate variables to use the saved
  final velocities (IREST=1, NTX=5) and letting it run twice as long
  (NTSLIM=2000):

	% sander -O \
		-i md1.in \
		-p stbtcap.top \
		-c md0.rst \
		-o md1.out \
		-x md1.crd \
		-r md1.rst


  Looking at the potential energy for both runs:

				Equilibration 2: Potential Energy Component

	


  The potential energy clearly has still not fully equilibrated. At this 
  point, it appears that the assumption that fast warming was ok could 
  have been wrong, and a gradual warming protocol or use of cartesian
  restraints in the initial warmup may have been in order (this more
  conservative approach is recommended in general). On the other hand,
  6-10 psec of equilibration is not an excessive requirement.

  In effect, we have performed a mild form of simulated annealing on the 
  initial structure, and in any case need to know how far the solute has 
  moved from the crystallographic starting position. For this analysis, 
  we turn to Carnal.

Carnal

  To get a rough measure of how far the structure has drifted from the
  crystal form during the equilibration, we perform a root-mean-square
  (RMS) comparison of the solute atoms that were allowed to move, using 
  the first set as a reference. (If this system were not held in place by
  the stationary atoms, we would want to get the best fit of each coordinate
  set to the first using Carnal's RMS FIT option.)

  Three types of RMS are measured: the entire moving part of the system;
  the moving biotin to get a rough idea of the changes that it experiences
  in its pocket; and the biotin using the FIT option to see how much of its 
  variation is due to internal (vs. orientational) changes.


	% carnal < carnal_rms.in > carnal_rms.out


  As might be expected, the average RMS decreases in each successive case:
  0.9, 0.7 and 0.2 Angstroms respectively. The time results:


				Equilibration: RMS Deviation from Initial Structure

	


				Equilibration: RMS Deviation of Backbone from Initial Structure

	


  Although the equilibration is not quite complete based on the potential 
  energy above, the overall RMS indicates that the structure is not drifting
  progressively (as virtually guaranteed by the frozen part). The overall
  RMS of ca. 1 Angstrom is well within the norm for simulations; beyond 
  2 Angstroms would be alarming. The backbone RMS is also reasonable:
  NMR-derived structures tend to differ from X-ray by about 0.8 Angstroms.
  (Different X-ray structures vary by about 0.4 Angstroms.)

  Turning to the biotin:


				Equilibration: RMS Deviation of Biotin from Initial Structure

	


				Equilibration: Internal RMS Deviation of Biotin

	


  The first graph indicates a quasi-periodic motion relative to the initial
  structure that may be interesting to explore once the trajectory has 
  equilibrated. The internal variations are faster and smaller, as one
  would expect. 

  Since the solute RMS is not drifting appreciably, a progressive
  change in water structure presumably is the cause of the progressive
  dropping of potential energy. This stands to reason, since we simply
  superimposed a sphere cut from a periodic system of pure water onto
  the solute and trimmed away any waters that overlapped, then warmed
  the local system rapidly to 300K, with nothing to prevent the sphere
  of waters from expanding. Therefore one would expect the drop in
  potential energy to correlate with the waters settling in toward the
  solute. Measuring the distance between the geometric centers of two
  groups of atoms, the biotin and the water oxygens, illustrates:

	% carnal < carnal_dist.in > carnal_dist.out

  resulting in md01.dist.

  Treating the trajectory as if it were at equilibrium for didactic purposes,
  we run carnal again to analyse hydrogen bonding between the ligand and its
  receptor:


	% carnal < carnal_hbond.in > carnal_hbond.out


  The summary data from the .out file shows that 6 hbonds persist throughout
  the trajectory ('grep 100 carnal_hbond.out'). They are:

    # 63 (SER  13 OG  )_(SER  13 HG  )..(BTN  479 O3  ) % occupied: 100.000000
    # 107 (TYR  29 OH  )_(TYR  29 HH  )..(BTN  479 O3  ) % occupied: 100.000000
    # 146 (ASN  35 N   )_(ASN  35 H   )..(BTN  479 O2  ) % occupied: 100.000000
    # 353 (SER  74 OG  )_(SER  74 HG  )..(BTN  479 O1  ) % occupied: 100.000000
    # 118 (BTN  479 N2  )_(BTN  479 HN2 )..(ASP  114 OD2 ) % occupied: 100.000000
    # 244 (BTN  479 N1  )_(BTN  479 HN1 )..(SER  31 OG  ) % occupied: 100.000000

  Thus all of the potential hbonding atoms of the biotin except for the sulfur
  are 100% engaged with the acceptor.  As one might expect, these hbonds are
  not weak; e.g. the first in detail:

    # 63 (SER  13 OG  )_(SER  13 HG  )..(BTN  479 O3  ) % occupied: 100.000000
        distance    avg       2.592  dev    0.113  max    2.971  min    2.410  N 150
        angle(deg)  avg       7.999  dev    4.300  max   21.185  min    0.890  N 150