The Frasch Lab          


    One project in the lab is using crystal structures to create computer simulations of the F1.  The computer simulations provide a tool that can help identify which residues are likely to be important in contributing to the function of the F1.  There are approximately 52,000 atoms in the F1 subunit making simulations extremely computationally intensive.       


This simulation is composed of 96 PDB files that were generated using 1BMF, 1H8E, 1E79, and 1E1Q, which are publicly available at the PDB website.  Another view of the simulation that shows specific residues involved in the escapement mechanism and the induction mechanism of the protein are available.  Simulations that include all 327,000 atoms involved in the process are currently underway.


 Comparison Between Adiabatic Mapping and Molecular Dynamics Calculations

Adiabatic Mapping

•         Potential Energy minimization is used to determine the atom positions for intermediate

        steps between one conformation and the next.

•         Need at least two existing conformations

–        Conformations must have identical atoms

•         Setup

–        Obtain/modify/create PDB files needed to describe protein conformations

–        Submit to Morph server

•         Number of frames between conformations

•         Provides quick results for gross protein motion

–        After submittal, receive results within 1-2 days

•         Atomic motion is not as accurate as full-blown Molecular Dynamics

–        How much off?

–        Not sure.  Don’t have any direct comparisons yet.


Molecular Dynamics

•         Computes trajectories by solving equations of motion using empirical force fields.

–        CHARMM

–        AMBER

•         Need one PDB file of protein placed in water box and charge neutralized

•         Need protein topology and parameter files

•         NAMD configuration file

–        Multiple parameters that tell NAMD about the protein

•         Setup

–        Create PDB file and PSF file

–        Submit to NAMD on cluster

•         Provides the most accurate in-silico results

•         Time step needed to accurately capture atom motion is 1 femtosecond (10-15 seconds)

•         Must be in the ms range to see any conformational changes in the catalytic site

–        1x10-3 s X 1x1015 fs/s x 1 step/fs = 1012 steps

–        Runs have required 1 week per million steps



Adiabatic Mapping

Two known F1-ATPase crystal structures (1E1Q and 1H8E) were modified and submitted to the Morph Server to obtain a series of files holding the atomic coordinates for each step of the morph as the protein changed from the first structure to the final structure.  The files were assembled within a TCL script to create a VMD movie showing the protein conformation changes occurring during 360° of g subunit rotation. Protein data file modifications needed before submittal to the Morph Server include removal of inhibitors used to obtain the crystal structures, addition of residues to create two conformations containing identical atoms, and matching of atom and residue names between the protein data files.

The current Yale morph server does not support the use of non-standard residues such as prosthetic groups, inhibitors, solvent molecules, and ions.  To complete the adiabatic mapping, the Morph server uses parameter files that describe the properties of standard amino acid residues.  Since the number of heteroatoms is in the tens of thousands and setting up parameter files for all possibilities is problematic, the authors of the Morph server have chosen not to include them.  Before any protein files are submitted to the Morph server, heteroatoms must be removed.  ATP, ADP, and Pi are heteroatoms and cannot be included in either of the F1-ATPase protein data files submitted to the Morph server.  For visualization purposes, ATP, ADP, and Pi atomic positions were estimated and added manually within VMD to each frame returned by the Morph server.

    Although atomic motion obtained from the Morph server may not be as accurate as full-scale molecular dynamics, the results “give a plausible representation for the pathway of the motion.  These movies can immediately give the viewer an idea of whether the motion is a rigid-body displacement or involves significant internal deformations…”  Since full-scale molecular dynamics requires extremely small timesteps (~10-15 seconds per step) to accurately capture molecular motions (Table 1), clustered parallel processors are needed to reach even picosecond (10-12 seconds) time scales after several weeks of continuous processing.  The Morph server provides results covering nanosecond to millisecond time scales after less than one day of processing time.  Because F1-ATPase requires milliseconds to hydrolyze ATP, the Morph server has provided insight into the molecular movements of its subunits much more quickly than molecular dynamics.

Table   1.  Molecular Timescales.

Distance and Forces between atoms in hydrogen-bond-forming residues

A trajectory movie of ATP hydrolysis in F1-ATPase created using the Morph server was analyzed within VMD to determine bond distances between residues in the g subunit and their counterpart residues in the a and b subunits.  We had identified these residues previously as important for rotation of the g subunit from mutagenesis studies (1,2).  Through TCL scripts similar to that shown in Figure 2, the center coordinates of each nitrogen or oxygen atom in the reactive portion of the important g subunit residue and the center coordinates of their bonding partner oxygen atoms in the reactive portion of a and b subunit residues was obtained for each residue-residue pair.  The TCL script steps through each frame of the simulation to obtain atomic center coordinates as the atoms move along their trajectories. 


Equation 1:                                   


Equation 2:                                                                            


Equation 3:                                                                                                               


Figure 3.  Atomic Distance and Coulombic Potential Equations.

    These center coordinates were imported into Excel to calculate the distance between bonding atoms and the associated force of electrostatic attraction, which follows Coulomb’s Law, as shown in Figure 3.  Equation 1 is the equation for the three dimensional distance between two points.  Equation 2 is Coulomb’s Law for bond energy, and Equation 3 is the equation for electrostatic force between two charges (or partial charges) and is derived from Equation 2 using a dielectric constant of 20 and lumping together all of the constants.  Finally, the x, y, and z components of each force were determined.  Within VMD, F1-ATPase is oriented so that the g subunit rotates around the x-axis.  This means that the important components of the forces acting on g are only in the y and z directions, and the final magnitude of each force is determined by summing the y and z vector components of the force.  Force vectors are shown as arrows.  Those vectors that act to contribute to rotation of g in the hydrolysis direction (counterclockwise when viewed from the foot of gamma) were defined to be positive (green arrow) when the atoms are approaching each other, and negative (red arrow) when they are moving apart (Figures 4 and 5).

Click on any of the following pictures to view a movie.


Figure  4.  Positive Force.                               Figure  5.  Negative Force.



The distance between potential hydrogen bonding atoms (mostly nitrogens and oxygens) was extracted using VMD and the corresponding electrostatic force was calculated within Excel for each residue pair at each step in the atomic trajectory.  Figure 6 shows the changes in distance as g rotates for g-Arginine 245 approaching and forming a salt bridge with b-Aspartate 316 (as shown in Figure 4).  Figure 7 shows the corresponding changes in electrostatic force calculated using Equation 3.  We hypothesize that the force contributing to rotation  increases to a maximum during the first 120 degrees as g-R254 approaches b-D316 and forms a salt bridge.  When the atoms involved in making the bond are at their closest the force is at a positive maximum.  As the gamma subunit continues to rotate, the force begins to act counter to rotation, thus becoming a maximum negative force.  We hypothesize that with  g-R254 and b-D316, it is at this point that b-D316 breaks the bond with g-R254 and forms a new salt bridge with a-R291, thereby minimizing the negative contribution of the  g-R254 and b-D316 interaction on rotation.

The sum of several residue interactions as a function of g rotation was created to estimate the net contribution of multiple residue interactions to generate rotational torque.  Figure 8 is a distance plot for some of these important residues.  Figure 9 is an individual force plot for these same residues.  Finally, Figure 10 shows the summation and total effect of these forces.

Figure 6.  Change in Distance between g-R254 and b-D316 versus g rotation angle.


Figure 7.  Change in Force between g-R254 and b-D316 versus g rotation angle.

Figure 8.  Distance plot of selected residues.


Figure 9.  Force plot of selected residues.

Figure  10. Summation of Forces for Residues shown in Figure 14 and 15. 

 Identifying potentially important residue interactions not found in crystal structures

The power of 3-D molecular modeling is evident when using the trajectory results from the Morph server.  An investigator may be able to infer important residue interactions by comparing crystal structures available from the Protein Data Bank.  However, the crystal structures are only snapshots of the molecule.  Having a simulation of the atomic motions as the conformation changes from one structure to the next provides insight into the potential interactions that may occur during the full motion of the protein.

For example, the residue a-Serine 335 does not show any obvious interactions with g in the crystal structures, but when evaluated within the simulation, it seems to be a very important residue that hydrogen bonds with g-Arginine 254.  An initial pass through the simulation has identified over eighty g residues that may hydrogen bond to (ab)3 residues. 


Although several residues important for ATP hydrolysis and g rotation have been identified, these do not completely explain g rotation.  Several more residues are thought to be involved in pulling g through its complete 360° rotation.  Using TCL scripts (Figure 6), g-a and g-b hydrogen bonding partners in each trajectory step can be identified.  Figure 11 shows as red the potential groups that may form hydrogen bonds and salt bridges between the g subunit and the (ab)3  ring as the groups come within 5 Angstroms of one another.

  Click on any of the following pictures to view a movie.

Figure   11.  Identifying Interacting Residues using VMD.

Several movies of F1Fo and F1-ATPase have been created using the trajectories determined by the Morph server.  These include:

·        ATP synthesis by F1Fo showing ADP and Pi moving into the bempty catalytic site, which disrupts the bonds between the g subunit and (ab)3, allowing the proton gradient across the c subunit of Fo to rotate both the c and g subunits.  As g rotates, changes in the conformation of (ab)3 lead to the conversion of ADP+Pi into ATP and a decrease in ATP affinity at another catalytic site leading to the release of ATP from that site (Figure 12)

·        ATP hydrolysis by F1 showing a closeup of ATP moving into the empty catalytic site, leading to disruption of g-(ab)3 bonds and conformational changes (Figure 13).

  Click on any of the following pictures to view a movie.

Figure   12.  F1Fo-ATPase Synthesizing ATP from ADP+Pi.


Figure   13.  F1-ATPase with ATP moving into empty catalytic site.



SwissProt- DeepView Protein Viewer (

Morph server (

The morph server interpolates between two protein conformations.  The user provides the conformations as PDB files and indicates the number of frames between the two conformations.  The process may take anywhere from five minutes to half an hour to complete, depending on the size of the structures and the load on the server. 

The interpolation is accomplished through adiabatic mapping, which minimizes the potential energy of the system at each intermediate step (frame) to determine the positions of all atoms for that step.


NAMD is “a parallel, object-oriented molecular dynamics code designed for high-performance simulation of large biomolecular systems.  NAMD uses the CHARMM force field and file formats compatible with both CHARMM and X-PLOR.  NAMD supports both periodic and non-periodic boundaries with efficient full electrostatics, multiple timestepping, constant pressure and temperature ensemble simulation methods.”  Depending on the size of the biomolecular system and the time scale of interest, results are obtained within a few days to several months or even years.


VMD is a 3D molecular visualization program that can be used to view, animate, and analyze molecular systems.  Any single protein file or series of protein files, whether obtained by morphing or molecular dynamics, can be viewed in VMD.  In addition to basic built-in viewing, zooming, and annotating capabilities, VMD can use custom TCL ( scripts to manipulate and extract atomic properties and trajectories, measure molecular geometry, and create movies showing molecular motion.  VMD can also be used interactively with NAMD to manipulate substrates in real time as the molecular dynamics simulation proceeds.

Microsoft Excel and Visual Basic:

Visual Basic scripts were written to prepare files for submission.  Microsoft Excel was used to calculate the hydrogen bond force for each residue-residue interaction, to calculate the components of the force, and to plot the results for each frame of the simulation.

Visual Basic was used within Excel to manipulate and reorganize data, insert and populate new worksheets for each residue-residue interaction, and was also used to automatically write TCL scripts that were then used within VMD.


Molecular Modeling - Molecular Dynamics using NAMD

Molecular Dynamics computes atomic trajectories by solving equations of motion for each atom using empirical force fields that describe the residue-dependent properties of each atom.  CHARMM and AMBER are two well-known force field formats.  As mentioned earlier, molecular dynamics requires femtosecond time steps to capture molecular motion.  To obtain results in a reasonable amount of time (weeks to months), the analysis must be run on a parallel cluster computing system.  All molecular dynamics simulations for this study were completed using NAMD on the ASU/TGEN Saguaro Cluster, which contains 1048 processors.  NAMD requires that the protein to be modeled be placed in solvent (water) with a neutral charge.  Figure 14 shows F1-ATPase within a water box, and Figure 15 shows ions added to the model.

Figure   14.  F1-ATPase in Water Box.

Figure  15.  F1-ATPase with Ions to Neutralize Charge.



uNeed one PDB file of protein placed in water box and charge neutralized

uNeed protein topology and parameter files

uNAMD configuration file

–Multiple parameters that tell NAMD about the protein


–Create PDB file and PSF file

–Submit to NAMD on cluster

uProvides the most accurate in-silico results

uTime step needed to accurately capture atomic motion is 1 femtosecond (10-15 seconds)

uMust be in the ms range to see any conformational changes in the catalytic site

–1x10-3 s X 1x1015 fs/s

x 1 step/fs = 1012 steps

–Runs have required 1 week per million steps



Echols, Milburn, and Gerstein (2003). Nucleic Acids Res. 31:478-82.

Krebs and Gerstein (2000). Nucleic Acids Res. 28:1665-75.

Gerstein and Krebs (1998). Nucleic Acids Res. 26:4280-90.

"The Yale Morph Server (".

"All images were made with VMD/NAMD/BioCoRE/JMV/other software support. VMD/NAMD/BioCoRE/JMV/ is developed with NIH support by the Theoretical and Computational Biophysics group at the Beckman Institute, University of Illinois at Urbana-Champaign."  (

   Yoshida, M., Muneyuki, E. and Hisabori, T. (2001) ATP synthase--a marvellous rotary engine of the cell. Nat Rev Mol Cell Biol, 2, 669-677.

    Junge, W., Lill, H. and Engelbrecht, S. (1997) ATP synthase: an electrochemical transducer with rotatory mechanics. Trends Biochem Sci, 22, 420-423.

    Gerstein,M. and Krebs,W. (1998) A Database of Macromolecular Movements. Nucleic Acids Res., 26, 4280.

    McCammon, J. A. & Harvey, S. C. (1987). Dynamics of Proteins and Nucleic Acids. Cambridge UP

   Eisenberg, D. & Kauzmann, W. (1969). The Structure and Properties of Water. Clarendon Press, Oxford.

    Humphrey, W., Dalke, A. and Schulten, K., (1996) VMD - Visual Molecular Dynamics J. Molec. Graphics, 14.1, 33-38.

    Caddigan, E., Cohen, J.,Gullingsrud, J., Stone, J., (2005) VMD User's Guide Version 1.8.3.  Theoretical Biophysics Group, University of Illinois and Beckman Institute, 405 N. Mathews, Urbana, IL 61801. (

    Echols, N., Milburn, D., and Gerstein, M. (2003). MolMovDB: Analysis and visualization of conformational change and structural flexibility. Nucleic Acids Res 31: 478-82.

    Krebs, WG and Gerstein, M. (2000). The morph server: a standardized system for analyzing and visualizing macromolecular motions in a database framework. Nucleic Acids Res 28: 1665-75.

 Kale, L., Skeel, R., Bhandarkar, M., Brunner, R., Gursoy, A., Krawetz, N., Phillips, J., Shinozaki, A., Varadarajan, K., and Schulten, K. NAMD2: Greater scalability for parallel molecular dynamics. (1999)  J. Comp. Phys., 151:283-312.

 Bhandarkar, M., Brunner, R., Chipot, C. , Dalke, A., Dixit, S., Grayson, P., Gullingsrud, J., Gursoy, A., Hardy, D., Humphrey, W., Hurwitz, D., Krawetz, N. Nelson, M., Phillips, J., Shinozaki, A., Zheng, G., Zhu, F. (1999) NAMD User's Guide, Version 2.5. Theoretical Biophysics Group, University of Illinois and Beckman Institute, 405 N. Mathews, Urbana, IL 61801.

  Per the NAMD license agreement, "NAMD was developed by the Theoretical and Computational Biophysics Group in the Beckman Institute for Advanced Science and Technology at the University of Illinois at Urbana-Champaign."

  Braig K, Menz RI, Montgomery MG, Leslie AG, Walker JE (2000) Structure of bovine mitochondrial F1-ATPase inhibited by Mg2+ADP and aluminium fluoride. Structure 8:567–573

  Menz RI, Walker JE, Leslie AG (2001b) Structure of bovine mitochondrial F1-ATPase with nucleotide bound to all three catalytic sites: implications for the mechanism of rotary catalysis. Cell 106:331–341

  Moss, D. (Course Organizer) (1998) The Principles of Protein Structure: Using the Internet '97 - '98, A Birkbeck College (University of London) Accredited Advanced Certificate Course. (


Last updated 2/20/2006

Frasch Lab

Arizona State University