Sophia Yancopoulos,  Mihaly Mezei, Roman Osman

Originally  presented as poster 210 February, 2001 at  Annual Biophysical Society Meeting  in Boston, MA

click for MMC



The recent crystallographic structures of the D96N mutant of bacteriorhodopsin (Luecke et al, 1999)in the ground (BR) and the excited (M) states reveal internal waters with proposed specific functionality. We have investigated the localization of the internal waters with grand canonical ensemble Monte Carlo simulations. The most probable positions of simulation sites were determined with the generic site method. The results of the simulations show 3 clusters of internal sites in each state and additional isolated sites in hydrophobic pockets. Ten of the 11 internal BR state waters are well matched by the simulations at an occupancy-weighted average distance of 0.69 ± 0.5 A. Three additional internal sites with total occupancy of 1.71 have been found for the BR state and 11 additional in the M state with total occupancy of 6.4 waters. Superimposing the two data sets yields 14 water pairs with an average distance of 0.82 ± 0.45 A. The simulation locates only one of three missing waters near Asn-96 in the M state.Supported by Training Grant DK07645


Water is an in integral component of proteins and plays an important role in their structure and function. High resolution structures determined by x-ray crystallography have recently become available for bacteriorhodopsin, a 7 TM protein, and its mutants (Luecke et al, 1999). The waters in these structures have been clearly observed and offer an excellent opportunity to test theoretical methods of solvation. Furthermore, simulations offer an opportunity to investigate the properties of internal waters and their role.


In order to explore the placement of water in the transmembrane domain of the 7 TM protein bacteriorhodopsin we have used a Grand Canonical Ensemble Metropolis Monte Carlo (MMC) cavity biased insertion method (Resat and Mezei, 1994). This method has already been successful in determining water positions in a number of crystal structures, for example,  small polydisaccharides (Resat & Mezei, 1994),  DNA (Resat & Mezei, 1996), and  complexes of inhibitors bound to HIV-1 protease (Marrone et al.1998). The cavity-biased formulation of MMC was chosen since it improves statistical sampling efficiency (Mezei 1987) and is therefore an appropriate method to investigate the solvation of the TM domain of BR.


In this approach the insertion of a molecule is attempted if a cavity of appropriate radius is found and accepted with a probability:

P' = min { 1,Pcav(N) exp[-(B+E(N+1)- E(N))/kT+B] / (N+1)}

where E(N) is the potential energy of the system of N particles and Pcav(N) is the probability of finding a cavity of a specific size.

The parameter B is related to the excess chemical potential as

µ' = kT (B - ln <N>)

Where <N> is the average over the number of particles.


The scheme is implemented by first adjusting the  B parameter  until the targeted number of molecules of solvent is approximately achieved outside a rectangular region containing the  molecule so that the density in this region external to the protein approaches the density of water (0.997 grams per cc).  Water molecules are then inserted or deleted by accepting or rejecting each attempt based on the Metropolis Monte Carlo criterion.

The tuned values of the B parameter were respectively: -1.2 and ?1.35 for the BR and M states.


The atomic coordinates of the BR (unphotolyzed) state, observed at 1.8 Angstrom resolution, and M (excited) state observed with a resolution of 2Å, from the Protein Data Bank (accession numbers 1C8RR and 1C8S respectively). The waters and lipids were removed from the structures.  Hydrogens were added to the structures with CHARMM. The N-and C-termini as well as two residues, which were discontinuous in the original structures, were capped with acetyl and N-methyl groups, respectively.

The BR and M states were run for  5 million and 6 million MMC steps, respectively, with the equilibrated B parameter. The coordinates of the atoms in the systems were stored every 10,000 steps during sampling.  During the simulations the solute residues were kept fixed at their experimental crystal structures and only the water molecules were allowed to move. This was in keeping with prior MMC simulations (Resat and Mezei 1996). Periodic boundary conditions were applied with the cell size as determined by simulaidês OPT. The temperature was 298 K. The water model chosen was the TIP3P (Jorgensen et al., 1983) with the solute potential type calculated in CHARMM.  The solute-solvent cutoff was used with the minimum-image convention. The solvent spherical cutoff was 10 Å.  The resultant average densities obtained in the exterior region were  0.99875 g/ml and 0.99743 g/ml respectively  in the BR and M states.


To establish the positions of water sites from the simulation,  we utilized the generic solvent (GS) site method developed by Mezei and Beveridge, 1984.  The configurations from the trajectory are superimposed and the generic sites are determined by the graph theoretical  method, (Berge, 1962), which minimizes the root mean square deviation between the GS and the individual waters assigned to it in each snapshot.  The GS are characterized by  mean locations for the oxygen positions, the mean square deviation of the position, and the occupancy of the site: the ratio of the number of water molecules assigned to the location divided by the total number of configurations.


                                                                                                         Table 1a: Internal BR State Generic Sites matching  Crystallographic Waters

Figure 1.(blue) BR state simulation waters  superimposed on
           (red) X-ray positions.


Two factors influence the occupancy of a site. One is the potential energy of the water at a particular insertion site. The other is the B factor, which controls the equilibration between the excess chemical potential and the average number of particles in the GC ensemble.

  Energies of Representative Sites in k-cal/mol


Conclusion: Waters in low occupancy sites have significantly higher energies.

                        Occupancy as a function of B:


Increasing the excess chemical potential (i.e., increasing B) progressively increases the  occupancy by overcoming the unfavorable energy.

Simulations find Additional Waters:

                                                                                                      Table 1b: Internal BR Generic Sites withoutmatching Crystallographic Waters:


FIGURE 2: A snapshot from the simulation. BR state Crystallographic water 405 is shown in yellow. One match Simulation water 236 is to the left.  An additional match, site  256 is on the right .  Another two  nearby simulation waters, 251 and 264,  are not seen in this view, but their properties are summarized   in the table and below.



Additional waters are placed in plausible positions.

MMC is successful at locating potentially real water sites that were not observed crystallographically.



Back to the Mezei Lab home page