Modeling the Bacterial Photosynthetic Reaction Center. 2. A

May 25, 1999 - Ab initio and other computational studies of bacterial reaction center cofactors are usually performed using the observed (low-resoluti...
1 downloads 11 Views 165KB Size
Download PDF
4906

J. Phys. Chem. B 1999, 103, 4906-4915

Modeling the Bacterial Photosynthetic Reaction Center. 2. A Combined Quantum Mechanical/Molecular Mechanical Study of the Structure of the Cofactors in the Reaction Centers of Purple Bacteria Michael C. Hutter,† Jason M. Hughes,† Jeffrey R. Reimers,*,† and Noel S. Hush†,‡ School of Chemistry and Department of Biochemistry, UniVersity of Sydney, NSW 2006, Australia ReceiVed: August 13, 1998; In Final Form: October 26, 1998

Ab initio and other computational studies of bacterial reaction center cofactors are usually performed using the observed (low-resolution) X-ray structures. Unfortunately, these geometries are necessarily approximate and this can have dramatic influences on calculated properties. For example, the calculated energies of the four bacteriochlorophylls in Rhodobacter sphaeroides vary by over 160 kcal mol-1. To overcome this problem, a combined quantum mechanical/molecular mechanical (QM/MM) method is employed to optimize the structure of the special pair and other cofactors in the photosynthetic reaction centers of Rhodopseudomonas Viridis and Rhodobacter sphaeroides, while a purely MM model is used to refine the structure of the remaining protein environment. Specifically, the QM/MM optimizations are performed using a semiempirical AM1based formalism. After relaxation, the energies of the bacteriochlorophylls differ by only typical conformer relative energies, ca. 5 kcal mol-1. Another example of improved cofactor properties is the PL-PM interaction energy which has been predicted to be strongly repulsiVe at the X-ray structure but here is shown to be realistically attractiVe after optimization. After optimization, the distortions in the geometries of the cofactor are seen to be controlled by protein-cofactor interactions, and the cofactors on the L-side are all seen to fit more snugly together within the protein environment than do their M-side counterparts. Also, the 2a-acetyl group of PM for Rb. sphaeroides, for which hydrogen bonding to the protein is restricted, is predicted to form a weakly bound sixth ligand to the magnesium of PL; this is consistent with, but not obvious from, the X-ray structure.

Introduction The initial steps of photosynthesis in which sunlight is converted into chemical energy involves a charge separation process from a chromophore to a primary electron acceptor. Within purple bacteria, photooxidation occurs in a membranebound protein which has been characterized by spectroscopic and crystallographic studies.1-9 The prosthetic groups of the reaction centers of these bacteria consist of four bacteriochlorophylls, two bacteriopheophytins, and two quinones, all organized in pseudo-C2 symmetry. This is illustrated in Figure 1. Charge transfer is found to occur very rapidly from the two central bacteriochlorophylls, which are referred to as the special pair, toward the quinones with excellent quantum yield.10 Despite the apparently symmetrical arrangement of the cofactors, the electrons are transported along only the L-branch of the reaction center, suggesting that there are interactions with the protein which favor this path. To account for this observed asymmetry, hydrogen bonding, electrostatic, and steric influences have been investigated.6-9,11-15 Analysis of ESR related spectroscopic studies has shown that this already affects the one-electron oxidation of the special pair leading to an unequal spin distribution on the M- and L-half.1,2,4 Up to now, a considerable number of computational investigations on the photosynthetic reaction center have been published employing various levels of theory.8,13,14,19-28 Owing to the large size of the surrounding protein, several groups have * To whom correspondence should be addressed. E-mail: reimers@ chem.usyd.edu.au. FAX-No.: +61 2 9351 3329. † School of Chemistry. ‡ Department of Biochemistry.

Figure 1. The prosthetic groups (cofactors) in the photosynthetic reaction center of Rh. Viridis.

used the finite difference Poisson-Boltzmann equation to treat its electrostatic influence.7,8,14,29 An explicit treatment of atomcentered charges was applied by Thompson and Schenter24 who also included polarizable dipoles in the protein region. Warshel, Parson, and co-workers30,31 used the LPDP (protein dipoles, Langevin dipoles) approach to calculate electrostatic energies of the cofactors and selected residues of the protein. For the investigation of spectroscopic properties, usually the INDO/S method32 is applied.13,19-23 Recent studies in our group employed INDO/S in multiconfiguration self-consistent-field (MCSCF) mode with multireference (singles) configuration interaction (MRCI-S) to describe the spectrum of the radical cation of the special pair.25,26,33

10.1021/jp9833808 CCC: $18.00 © 1999 American Chemical Society Published on Web 05/25/1999

Bacterial Photosynthetic Reaction Center While the above studies apply very sophisticated treatments of properties, such as the electronic structure of the cofactors and/or the protein-cofactor interactions, typically the geometries used for the cofactors were simply taken from the X-ray data (adding the required hydrogen atoms by heuristic principles). In contrast to crystallographic structures of small organometallic compounds, however, the corresponding atomic positions in the reaction center protein are less well-resolved. As computed properties are very sensitive to geometry, severe errors have been introduced. A typical example of this is the relative ab initio SCF energies of the bacteriochlorophylls PL and PM: evaluated at the X-ray geometry in Rh. Viridis,29 these differ by a staggering 120 kcal mol-1, a value in error by ca. 2 orders of magnitude. Clearly, any computed result obtained using such crude nuclear coordinates must be heavily scrutinized for accuracy. To date, very few calculations have been reported that optimize or refine the molecular structure of the bacteriochlorophylls or bacteriopheophytins by the means of quantum mechanics with respect to the surrounding protein. There have been various molecular dynamic studies using standard force fields which have considered relaxation and entropy of the entire protein.7,8,11,28,31 After this work had been submitted, Siegbahn and co-workers18 reported optimizations of the special pair and other cofactors by using a “piecemeal” approach using the hybrid density-functional type methods with essentially double-ζ quality basis sets in the gas phase. That is, the strategy they employed involved breaking the reaction center into small subunits to form medium sized molecules which ensured computational tractability. The surrounding protein was treated as a dielectric medium and the cavities around the solute molecule were determined by isodensity surfaces. Essentially they were concerned with modeling the electron transfer processes in the reaction centers of Rhodobacter Sphaeroides and of photosystem II where good agreement (within 2 kcal mol-1) with the experimental separation energetics was found. They also highlight the role of the local protein environment in respect to tuning the properties of the various chromophores. Also, Hasegawa and co-workers16,17 used the symmetry adapted cluster-configuration interaction (SAC-CI) methodology on the reaction center of Rhodopseudomonas Viridis to study the mechanism and unidirectionality of the electron transfer from the photoexcited special pair to the bacteriopheophytin. In addition, they calculated the theoretical SAC-CI spectrum for the ground and some excited states of the photosynthetic reaction centers. Although the protein, waters, and other cofactors where dealt with using the point-charge electrostatic model, the geometries of the special pair and cofactors, however, were based on the X-ray structure. While quantum-mechanical methods offer significant advantages over force-field calculations, particularly in relation to the structure and properties of the cofactors, they are not currently feasible for studies of the entire protein. In this study we model the in situ structure of the cofactors of Rh. Viridis and Rh. sphaeroides. Two aspects are involved: the generation of protein structures and the optimization of the cofacter geometries within these structures. Initial protein structures are obtained starting with the X-ray structures, adding hydrogens and correcting for known crystallization effects. These are then relaxed using the Amber molecular mechanics force field. Accurate cofactor structures are then obtained using a combined quantum mechanical/molecular mechanical (QM/ MM) technique which combines the AM1 and UFF49 methodologies.

J. Phys. Chem. B, Vol. 103, No. 23, 1999 4907

Figure 2. Chromophores in the investigated bacterial reaction centers.

The Initial Protein Environment The photosynthetic reaction center of the purple bacteria studied here, Rh. Viridis and Rb. sphaeroides, consists of several protein subunits from which the L and M amino acid sequences form trans membrane helices while the cytoplasmic H chain shows far less secondary structuring.5,6,9 The cofactors of both bacteria differ only slightly. As shown in Figure 2, the chromophores in Rh. Viridis are bacteriochlorophyll b and bacteriopheophytin b while in Rb. sphaeroides the corresponding a types are present. Rb. sphaeroides contains two ubiquinones while one of these, QA, is a menaquinone in Rh. Viridis. A nonheme iron cation is present in both bacteria. Each contain a carotenoid on the M-side; this is spheroidene in the wild type strain of Rb. sphaeroides but is dihydroneurosporene in Rh. viridis. The high sequence homology for the protein subunits and especially the very similar overall secondary structure of both reaction centers highlights the functional similarity of both reaction centers.6,9 Starting with the X-ray coordinates of Rhodopseudomonas Viridis (1PRC)5 and Rhodobacter sphaeroides (1PCR)6 as given in the Protein Data Bank at the Brookhaven National Laboratory,53 hydrogen atoms and atom centered charges were added to the amino acid residues using template values according to the AMBER force field.54 For this purpose we used the HYPERCHEM program,55 adding hydrogen atoms also to the cofactors. Unfortunately, various issues arise concerning how the observed X-ray structures relate to the actual situation in ViVo, and these are addressed in turn below. As described later, all ionizable residues in the hydrophilic region of the protein are assumed to be charged (no counterions are added), while all those in the hydrophobic region are assumed to be charge neutral. Missing from the X-ray data are various terminal amino acid residues. As these are far removed from the sites of interest within the protein, no attempt is made to replace these missing groups, and the chain ends are also capped with H atoms (according to valency). Conversely, present in the X-ray data are some lauryldimethylamine-oxide molecules which were introduced during the extraction process. To make the protein model more representative of the situation in ViVo, these molecules are removed. Because the C chain of Rh. Viridis (which contains four heme groups) is not found in the corresponding X-ray data of Rb. sphaeroides,6,9 it is omitted. This provides a significant reduction

4908 J. Phys. Chem. B, Vol. 103, No. 23, 1999 in the required computational effort and facilitates comparison of EvdW for the two bacteria. Except for possible electric field effects and a proposed role in re-reduction of the photooxidized special pair, no significant effects of the C chain on the structure and function of the remainder of the reaction center are expected.57,58 Various numbers of water molecules are found in the X-ray structures. For Rb. sphaeroides, all 160 crystallographic water molecules are included. However, for Rh. Viridis, many of the 201 crystallographic waters are not in close contact with the amino acid chains; we removed these, leaving a total of 150 waters (an alternative treatment39 would be to add additional waters to complete the structure). It is thought that the water molecule orientation may be very important in the electronproton transfer mechanism as they may give rise to several possible water channels (proton paths) connecting the secondary ubiquinone to the aqueous phase.59 Due to the fact that X-ray structures do not determine the orientation of the water molecules in the protein, the HYPERCHEM program55 was used to orient the waters with respect to the protein environment (in this initial step, the heavy-atoms’ coordinates were held fixed and the temperature was at 500 K). For all water atoms, TIP3P charges60 were assigned (0.417 e for hydrogen and -0.834 e for oxygen). Using the geometry optimiser in HYPERCHEM55 and employing the steepest descent method in conjunction with the Amber molecular mechanics force field, the waters were reoriented until the root-mean-square gradient reached a value that was less than 0.1 kcal Å-1 mol-1. A nonheme iron atom is found in both bacteria and is usually assumed5,6 to be in a formal oxidation state of II, although it is sometimes listed as Fe(III).7 Also it is known that the iron does not change its formal oxidation state in the electron transfer process and even substituting the iron atom with that of zinc (whose redox potential is out of range) the electron transfer process will still proceed. Thus, the nonheme iron atom was assigned a formal oxidation state of II for all calculations. It is ligated by four histidine (M219, M266, L190, and L230) and one glutamate (M234) residue in approximately octahedral coordination. The atomic charges for the sulfate (SO4 617 through to SO4 623 inclusive) and phosphate (PO4 800) groups found in the X-ray structures of Rhodopseudomonas Viridis and Rhodobacter sphaeroides, respectively, were determined by fitting the molecular electrostatic potential (ESP).56 The derived atomic charges were calculated at the RHFESP/6-31 G* level of theory for the sulfate (S ) +1.407 e, O ) -0.852 e) and phosphate (P ) +1.512 e, O ) -1.128 e) groups. The close vicinity of the carboxyl side chain of GluL104 to the 9-keto group of ring V of BPhL found in the X-ray structures in both bacteria suggests the presence of a hydrogen bond.5,6,9 For this to occur, the residue must be protonated making it electrically neutral, and a hydrogen atom was added. Charges for this protonated residue were obtained by using the RHF/631G* level of theory and the calculated ESP charges for butyric acid were assigned to be 0.505 e for the new proton and -0.770 e for the oxygen of the hydroxyl group. Applying the AMBER partial charges of the GLN motive for the remaining atoms of this amino acid yields -0.384 e for the keto oxygen of the carboxyl group, the residue becoming net neutral. For most residues (especially histidine) capable of undergoing acid-base reactions, the actual extent of protonation is uncertain. The reaction center protein contains both hydrophobic and hydrophilic regions, with the ionizable residues located almost

Hutter et al. entirely in the hydrophilic part; located therein, the pKa of the amino acid side chains at physiological pH is such that all ionizable residues are most likely present in their charged form, and this is assumed in this study. Nevertheless, other studies of the reaction center have used different assumptions, and there is no generally agreed recipe. Gunner et al.,14 for example, used a protein model for Rh. Viridis which included the C chain and had a total charge of -6 e, but in this study all of the histidine residues and the sulfate groups were taken to be “uncharged”. Other studies,24,30 however, have altered the atomic charges of the ionizable amino acids to make each individual residue electrically neutral. Because of the present uncertainties, it was decided better not to force exact charge neutrality. The net charge for Rh. Viridis and Rb. sphaeroides are -1 e and +2 e, respectively, charges which are very small for proteins of this size. In the X-ray structure of Rh. Viridis, seven sulfate groups are reported, and it is likely that these entered during the crystallization process. For Rh. Viridis only the sulfate 618 is retained since it is ligated by two positively charged residues (ArgM265 and HisM143). On the other hand, Rb. sphaeroides contains a (noncrystallographic) phosphate ion (800). It arises due to the presence of eight additional positively charged residues on the (hydrophilic) H chain of Rb. sphaeroides. Also, for this bacterium, one negatively charged residue (GluH258) is not given in the crystal structure. QM/MM Methodology The combined QM/MM method used in this study is described in full detail in ref 34 and is referred to hereafter also as PCM (point charge model). The corresponding option is implemented in the VAMP program package35 which has been used for the majority of the calculations. It has been successfully applied to perform docking,36,37 inhibitor-strength ranking,37 and the calculation of heats of adsorption in zeolites.34,70 Within the formalism of the PCM, the molecular systems are partitioned into quantum mechanical (QM) and molecular mechanical (MM) parts which are completely separated from each other. This separation is straightforward if no bond breakage is involved (i.e., molecules are included entirely only in one part), but this is not the most common situation and a variety of separation schemes have been introduced.36-45 In this study we consider only the bacteriochlorophyll, bacteriopheophytin, carotenoid, and other cofactor molecules in the QM part. Separation of the QM and MM parts is thus straightforward except for the bacteriochlorophylls which are attached to the protein backbone via histidine ligands (L173, M200, L153, and M180 of Rh. Viridis and L173, M202, L153, and M182 of Rb. sphaeroides for PL, PM, BChlL, and BChlM, respectively). Our approach is to chop the imidazole side chains from the protein backbone, including them with their bacteriochlorophylls in the MM part. This is achieved by deleting the histidine β CH2 units and terminating the γ carbons with H atoms. The R C atoms remain in the MM part, and their atomic charge is renormalized so that the net charge of the MM part is not modified. Within the PCM formalism, the larger molecular part (the protein environment) is described in a force-field-like fashion while the smaller part (the cofactor(s) being considered) is subject to semiempirical molecular orbital treatment based on AM146 theory. For the bacteriochlorophylls, the AM1 magnesium parameters published in the proceeding paper of this series47 are used. The total energy of the system is given by the sum of the QM part, the MM part, when necessary a term arising

Bacterial Photosynthetic Reaction Center

J. Phys. Chem. B, Vol. 103, No. 23, 1999 4909

from the QM-MM separation, and the interaction energy arising from the MM part toward the QM part:

k

∑γ 2 |δrγ|2

(2)

where γ indexes the Cγ atoms on each partitioned histidine and δr is the displacement of that atom away from its initial location. A value of 10000 eV Å-2 is used for the force constant, k, this value being sufficiently large to prevent the Cγ atoms moving more than ca. 0.003 Å from their initial positions. The QM/MM interaction energy consists of two terms including van der Waals and electrostatic contributions:

EQMrMM ) EvdW + Eestat

(3)

The van der Waals term contains a Lennard-Jones (12-6) potential which is similar to the expression used in force fields49 and other QM/MM studies.24,45,50

∑A ∑k

( [ ] [ ])

DAk -2

RAk

6

+

rAk

RAk

12

(4)

rAk

where A indexes the QM atoms, k indexes the MM atoms, DAk is the well depth for a specific pair of atoms defined as the geometric mean of the single atomic van der Waals energies DA and Dk, RAk is the equilibrium distance between these atoms, and rAk is the actual separation. The corresponding parameters of the UFF force field were used in this implementation.49 While the van der Waals interactions are subject to extensive parameterization, the electronic interaction energy which leads to a polarization of the QM part by the MM part is given by simple electrostatics which can be divided into core and electronic terms:

Eestat ) Ecore + Eelec ) ZAqk

(

qk

∑A ∑k 4π r ∑A ∑k ∑µ ∑ν 4π Pµν‚(µν|sksk) 0

Ak

0

∑k 4π (µν|sksk)

(6)

(1)

In this, EQM and EMM are the internal purely quantummechanical and molecular-mechanical energies, respectively. Here, a very simple form of the QM/MM separation energy (readily implemented in VAMP) is used, namely, the constraint48

EvdW )

qk

0

Etot ) EQM + EMM + ESEP + EQMrMM

ESEP )

Fµν ) Foµν -

)

(5)

where  is the effective dielectric constant of the medium, ZA are the QM atomic core charges, qk are the MM atomic charges, Pµν are elements of the QM density matrix for atomic orbitals µ and ν on atom A, and (µν|sksk) are two-center one-electron AM1-type integrals describing the interaction between µ and ν (on atom A) and an s-type function of zero expansion.34,51 Note that the form of these two-center integrals corresponds to the multipole expansions sets of Dewar and Thiel51 describing the interaction between a point charge and a multipole of arbitrary order, but here we include only contributions of up to second order. In particular, (µµ|sksk) ) 1/rAk and so the leading term in this equation is simply34 the Coulomb interaction of the MM charges with the net QM atomic charges ZA - ∑Pµµ. The calculation of the new Fock matrix is similar to other QM/MM studies which involve semiempirical Hamiltonians40,50 with the elements of the perturbed Fock matrix given simply by

where Foµν is the unperturbed QM-only contribution. Voityuk et al. have suggested that an additional parameter be introduced into these integrals, fitted so that calculated gas-phase hydrogenbond energies agree with observed data for the ion-molecule complexes.50 We do not do this and instead the two-center integrals (µν|sksk) are simply those as used by Wang and Ford in a SCRF study.52 During the QM/MM optimizations we keep the protein frozen and hence the energy of the MM part is not affected by structural changes in the QM part and does not require evaluation. The assumption of a rigid MM part is justified when environments which show only reasonably small nuclear movement upon interaction with the QM part. The adsorption of small organic molecules into a zeolite cage is a suitable example. However, this approach can also be employed for the description of enzyme-ligand complexes involving nonallosteric proteins.34 Compared to other methods (eg. refs 13, 14, 16-18, 24, and 29-31) used for modeling the protein-cofactor interactions, the electrostatic interactions are treated very simplistically in our PCM model. While the protein consists of hydrophobic regions surrounding the cofactors and distant hydrophilic regions, our model embodies only one dielectric constant . Specific polarization effects are not included, as are effects arising from boundary regions not represented in the X-ray structure. However, these factors should not significantly effect the geometrical structure of the cofactors. For the dielectric constant, a value of  ) 1 applied to the interaction of atoms k and A in eq 5 effectively assumes that all material located between these atoms in the MM part is unpolarizable. This is, however, only appropriate for the interactions between the cofactors and their nearest neighbor atoms in the MM part or if electronic and geometric polarization is explicitly included in the calculation. If polarization is not explicitly included, then all other charges in the protein are shielded by the inner molecules and a value of  > 1 is appropriate. Values between 2 and 4 or even distance-dependent screening functions of the dielectric constant have been proposed to account for electronic effects involved such as polarizability and nuclear relaxation14,39,61,62 and these have been employed in other studies for which long-range dielectric properties are of central importance.14,29 If the entire protein, membrane, and aqueous structures are known, then the nuclear relaxation contribution can be treated explicitly and the appropriate value of the dielectric constant becomes the high-frequency component ∞. Difficulties arise as not all of the relevant structure is seen by the X-ray experiments30,31 and hence  is also used to implicitly account for the missing pieces. These missing pieces include both hydrophobic parts, for which low dielectric constants are appropriate, and hydrophilic parts, for which they are not. Again, for the hydrophobic regions, however, (those surrounding the cofactors), the high-frequency dielectric constant ∞ is most appropriate. This is approximately equal to the square of the index of refraction and is about 2 for most organic liquids and water. The problem of assigning a suitable constant for proteins has also been reviewed by Harvey.63 For the described PCM formalism, it was shown that a dielectric constant of 2 yields the lowest deviations in predicting the length of hydrogen bonds compared to crystallographic data.34 Hence this value was employed throughout this study.

4910 J. Phys. Chem. B, Vol. 103, No. 23, 1999

Hutter et al.

TABLE 1: AM1 Energies, in kcal mol-1, of the Cofactors at the X-ray Structure (with Added Hydrogens) and at the PCM-Optimized geometries, as Well as the Energy Change on Optimization Rh. Viridis

Rb. sphaeroides

cofactor

X-ray

PCM

change

X-ray

PCM

change

PM PL BChlM BChlL BPhM BPhL QA QB

65 -50 115 -36 18 58 47 -83

-219.1 -215.3 -217.3 -218.5 -150.6 -151.2 -53.3 -107.3

-284 -165 -332 -183 -169 -209 -100 -24

-115 -56 -92 -58 -37 -42 77 25

-246.5 -243.7 -240.6 -245.5 -174.4 -175.3 -131.0 -134.4

-132 -187 -148 -187 -137 -133 -208 -159

Geometry Optimization Procedure The calculations were performed in essentially four major steps. Initially, as previously described, the water molecules were oriented using the molecular mechanical geometry optimization routine as found in the HYPERCHEM program.55 During this, the protein and photosynthetic reaction centers and its associated cofactors remained frozen. Next, the cofactors were optimized in a preliminary QM/MM step inside a fixed protein environment. Next, the entire protein environment including the waters were optimized while the bacteriochlorophylls (including the histidine ligand attached to the Mg centers) bacteriophyiophtins, cofactors, and carotenoid remained frozen. This step is deemed particularly important since the initial structure of the protein was sourced from X-ray crystallographic studies and hence any molecular mechanics force field or quantum mechanics methodology would undoubtedly find some unrealistic local geometries. It is also noted in the previous study by Siegbahn and co-workers18 that the protein environment is crucial in fine tuning calculated molecular properties of the photosynthetic reaction center. Here, HYPERCHEM55 was again used within the molecular mechanics geometry optimization framework. It was found that the most efficient scheme was to subject the protein to 1000 cycles of the steepest descent minimization followed by using the conjugate gradient method until the root-mean-square gradient was less than 0.1 kcal Å-1 mol-1 (the resultant structure was checked to ensure that no abnormally short bond lengths or bond angles were found). In the final step, the structural optimization of the prosthetic groups was completed in sections using the remaining cofactors and the now optimized protein environment in the described QM/MM method. The sequence of optimization was QB, QA, BPhM, BPhL, BChlM, BChlL, carotenoid, then special pair (P). After optimization is completed for a particular group, its new atomic positions and AM1 derived charges are used to update the MM part for the subsequent calculations. In total, two complete QM/MM optimization cycles were performed, but the second cycle produced insignificant modification. Results The optimized structure of the cofactors in the protein environment are given in Supplementary Data as ASCII “.hin” files containing hydrogens (rhvd.hin and rbsp.hin) for use with HYPERCHEM.55 Detailed comparisons between the crystal structures of both reaction centers are given in refs 6 and 9; here we focus on structural and energetic changes resulting from the geometry optimization. Calculated energies for the cofactors both before and after optimization are shown in Table 1 while geometrical properties relating to individual interactions, particularly hydrogen bonding, are listed in Tables 2 and 3.

TABLE 2: Key Intermolecular Interaction Lengths Involving Cofactors of Rh. Wiridis in angstroms cofactor site PM

PL

PM BChlM BChlL BPhM BPhL QA QB

other site Mg 2a-keto-O 10a-keto-O 2a-acetyl Mg 2a-keto-O 9-keto-O 2a-acetyl Mg Mg 9-keto-O Mg 9-keto-O 10a-ester-O 10a-ester-O 9-keto-O proximal keto-O distal keto-O proximal keto-O proximal methoxy-O distal keto-O distal keto-O distal methoxy-O distal methoxy-O

HisM200-NE2 TyrM195-OH SerM203-HG CH2 H - PL Mg HisL173-NE2 HisL168-HE2 TyrL248-HG1 CH2 H - PM Mg PL Mg HisM180-NE2 Water304-H HisL153-NE2 Water302-H TrpM127-HE1 TrpL100-HE1 GluL104-H HisM217-HD1 AlaM258-H HisL190-HD1 HisL190-HD1 SerL223-HG GlyL225-H GlyL225-H AlaL226-H

1PRCa PCM 1.890 2.115 3.287 2.367 2.225 1.780 2.612 3.087 7.579 2.091 2.038 2.166 1.232 1.742 2.090 1.757 2.134 2.078 1.723 2.397 3.046 1.749 3.256 3.426

2.144 1.967 2.060 2.494 2.117 1.948 2.018 3.068 8.312 2.100 2.130 2.103 2.232 2.117 2.825 2.064 2.185 2.081 2.417 2.163 3.520 2.233 3.002 3.772

a

Name as given in the Protein Data Bank at the Brookhaven National Laboratory. (See refs 5 and 6.) Hydrogens were added as described in the text.

TABLE 3: Key Intermolecular Interaction Lengths Involving Cofactors of Rb. sphaeroides, in angstroms cofactor site PM

PL

PM BChlM BChlL BPhM BPhL QA QB

other site Mg 2a-keto-O 2a-keto-O 2a-acetyl-CH2 H 9-keto-O 10a-keto-O Mg 2a-keto-O 9-keto-O 2a-acetyl-CH2 H Mg Mg 9-keto-O Mg 9-keto-O 10a-ester-O 10a-ester-O 9-keto-O 9-keto-O proximal keto-O proximal keto-O distal keto-O distal keto-O distal keto-O distal methoxy-O proximal methoxy-O

HisM202-NE2 TyrM210-OH PL-Mg PL Mg IleM284-1HD1 SerM205-HG HisL173-NE2 HisL168-HE2 MetL248-2HE PM Mg PL Mg HisM 182-NE2 Water62-H HisL153-NE2 Water59-H TrpM129-HE1 TrpL 100-HE1 TrpM252-HH2 GluL104-H HisM219-HD1 ThrM222-HG1 AlaM260-H IleL224-H TyrL222-HH IleL224-H Water65-H

1PCRa PCM 2.152 3.374 3.374 4.035 2.496 3.439 2.257 2.684 2.262 2.106 8.008 2.216 1.947 2.308 3.165 1.927 1.989 2.427 1.800 2.188 2.960 1.769 1.980 6.260 2.167 2.852

2.114 3.814 2.349 4.562 2.450 3.929 2.174 2.023 2.693 2.585 7.738 2.099 2.249 2.109 5.405 2.345 2.464 2.458 2.061 2.203 3.837 2.132 2.281 6.971 2.238 6.189

a Name as given in the Protein Data Bank at the Brookhaven National Laboratory. (See refs 5 and 6.) Hydrogens were added as described in the text.

Energetics. As shown in Table 1, the relative energies of similar cofactors at the X-ray geometry varies by up to 165 kcal mol-1, while at the relaxed geometries this variation is dramatically reduced to less than 5 kcal mol-1, of the order of the energies associated with conformational variations in the cofactors. The relaxation energy of the macrocycles varies from -132 to -332 kcal mol-1. It is expected that analysis of cofactor properties evaluated at the PCM geometries will provide

Bacterial Photosynthetic Reaction Center

J. Phys. Chem. B, Vol. 103, No. 23, 1999 4911

TABLE 4: Protein-Cofactor Interaction Energiesa Rh. Viridis

Rb. sphaeroides

cofactor no.b

EvdW

Eestat EQM-MM no.b

PL PM BChlL BChlM BPhL BPhM QA QB

-108 -106 -89 -82 -117 -100 -91 -35

-25 -28 -13 -18 -10 -11 -8 -9

147 147 147 147 139 139 112 36

-133 -134 -102 -100 -127 -111 -99 -44

EvdW

Eestat EQM-MM

149 -103 -14 149 -107 -17 149 -87 -13 149 -77 -24 141 -120 -14 141 -93 -10 114 -79 -12 114 -65 -8

-117 -124 -100 -101 -134 -103 -91 -73

a Energies are in kcal mol-1. b no. number of atoms included in cofactor. All bacteriochlorophylls have an attached axial imidazole included within the cofactor, and all cofactors other than the selected one are treated as part of the protein.

Figure 3. Overlay of the optimized cofactors (gray) and the X-ray data (black) for Rh. Viridis.

significantly more reliable results than those obtained at the original X-ray geometries. Rhodopseudomonas Wiridis Geometry. An overlay of the optimized cofactors (gray) and the X-ray data (black) is shown in Figure 3. The conformational very flexible side chains show the largest structural deviations, especially for the phytyl tails of BPh on the M side and the carotenoid. The arrangement as well as the individual conformations of the phytyl groups also differ significantly in the crystallographic structures of both bacteria.6 These differences arise as the Rh. Viridis phytyl tail resides in a pocket which is accessible to the membrane surrounding the protein. In the crystal structure this tail is not resolved; only template atomic coordinates are listed with the atoms being assigned either large temperature factors in excess of 60 Å2 or zero occupancy. Our QM/MM-optimized structure should thus be considered as only one of a variety of possible conformations for the chain, the actual conformation(s) being determined by membrane interactions. With the exception of this chain and the carotenoid as well as the ends of the phytyl chains in Rp. sphaeroides, for which similar arguments apply, all other conformations are closely conserved during the optimization leading to an overall heavy-atom root-mean-square (RMS) deviation of just 1.02 Å, 39% of the stated resolution of the X-ray structure. Compared to X-ray data, the hydrogen bonds on the PL-half are slightly weakened so that the interaction with TyrM195 on the M side is much stronger within the special pair. The intermagnesium spacing is increased by about 0.7 Å (Table 2), while the overlapping of ring I is unchanged (see below). The conformations of the 2a-acetyl groups are closely preserved except for BChlL for which it is rotated by ca. 30°. Hydrogen bonds involving ester keto oxygens are found to be less preferred than those involving carbonyl oxygens, particularly for the 10a ester groups of the bacteriopheophytins where these interactions

Figure 4. Overlay of the optimized cofactors (gray) and the X-ray data (black) for Rb. sphaeroides.

are weakened. The resulting distance (2.1 Å) between the 9-ketooxygen and the hydrogen atom of GluL104 supports the assumption that the carboxyl group of this residue is in fact protonated.5,6,9 The importance of orientating the waters found in the crystallographic structure is highlighted by the interaction of the 9-keto-O (located on the BChlL) with the hydrogen of water 302. For this the starting bond length is an unrealistic value of 1.2 Å, but after orientation, the bond length became realistic, 2.2 Å. For the menaquinone (QA), only minor variations of the hydrogen bonds were found while the ubiquinone (QB) is shifted away from the nonheme iron, averaging its interactions with the amino acids. The calculated axial magnesium-nitrogen bond lengths for the bacteriochlorophylls fall within a narrow range from 2.100 to 2.144 Å, while the corresponding X-ray distances vary from 1.9 to 2.3 Å. We investigated this bond type earlier using density functional and semiempirical calculations47 for a model magnesium porphyrin molecule with an axial imidazole ligand. Both the AM1 and B3LYP/3-21G methods predicted a Mg-N bond length of 2.14 Å, with the magnesium being about 0.29 Å above the ring plane. We thus expect that the calculated results for the reaction center are quite reliable: they suggest that the axial histidine residues are tightly bound to the bacteriochlorophylls. The porphyrin rings of the accessory bacteriochlorophylls (BChlM and BChlL) are nearly planar while the somewhat twisted conformations remains for those forming the special pair, most pronounced in the M-half in accordance with the X-ray structure. Displacement of the magnesium atom out of the porphyrin plane is found in similar five-coordinated complexes and is largest for BChlM and BChlL; for the special pair, the magnesium atoms are almost in plane with the nitrogen atoms, again in agreement with crystallographic data. Rhodobacter sphaeroides Geometry. The results of the cofactor optimization confirm a number of interactions discussed by Ermler et al. from their crystal structure6 and resemble the results for Rh. Viridis with an overall RMS deviation from the starting X-ray structure of just 0.90 Å; both structures are shown in Figure 4 while calculated structural properties are given in Table 3. The 9-keto oxygen of BPhL behaves analogously to that of Rh. Viridis, willingly adopting a hydrogen bond of length 2.0 Å to GluL104 and supporting the presupposition that the carboxyl group is not deprotonated. Some variation is seen between the calculated and observed intermagnesium separation within the special pair. The optimized values are of 7.7 Å and 8.3 Å for Rb. sphaeroides and Rh. Viridis, respectively, while the X-ray structures report 8.0 Å and 7.6 Å, respectively. Little variation is seen in the calculated axial magnesium-nitrogen bond lengths; however, these range from 2.099 to 2.174 Å. In contrast to Rh. Viridis, the hydrogen bonds involving the ester

4912 J. Phys. Chem. B, Vol. 103, No. 23, 1999

Hutter et al.

keto oxygens of the bacteriopheophytins in Rb. sphaeroides remain more significant, and the binding positions of the quinones are almost unaffected. One significant difference between the structures of both bacteria does emerge. This is consistent with, but not obvious from, the original X-ray data. For Rb. sphaeroides, a strong hydrogen bond is evident between the 2a-acetyl oxygen of PL and the hydrogen atom at HisL168, with the hydrogen-bond length reduced to 2.0 Å, and similar configurations for this group are also found for both PL and PM in Rh. Viridis. However, for the corresponding 2a-acetyl group of PM, hydrogen bonding to the protein would need to be via TyrM210 but this group is actually quite distant, see Table 3. Instead, the PCM calculations predict that the 2a-acetyl group of PM undergoes a conformational change which brings the keto oxygen to 2.4 Å away from the magnesium of PL, essentially forming a sixth ligand to the metal. Alternatively, the magnesium of PM and both centers in Rh. Viridis actually interacts weakly with the CH3 part of the 2a-acetyl groups, resulting in short Mg-H separations of ca. 2.4 Å, see Tables 2 and 3. SensitiVity to Long-Range Electrostatic Interactions. A central assumption used in the construction of our PCM model is that long-range electrostatic interactions do not significantly effect the optimized cofactor structures. We tested this hypothesis by using a very different model for the long-range electrostatics, one in which all histidine residues are in their neutral, unprotonated form; this model is more similar to the protein models used others, e.g., see refs 14, 24, and 30. We then reoptimized the structure for the cofactor which lies closest to the histidines, BChlM from Rh. Viridis. The result was a RMS change in coordinate of 0.08 Å and a RMS change in the internal bond lengths of just 0.0033 Å. Hence we see that the cofactor structures are insensitive to long-range electrostatic effects. Interpretation The Energetics of Dimer Formation. An important question concerning the reaction center is the nature of the PL-PM interaction, and the forces which hold the dimer together. If these are known, then it is, in principle, possible to design mutant reaction centers in which this interaction is controlled and hence the function of the reaction center modulated. Sakuma et al.29 have shown that, at the ab initio SCF level, an energy of ca. 21 kcal mol-1 must be supplied, at constant Rh. Viridis PL and PM intramolecular geometry, in order to push the two halves together. This compression energy is supposedly provided by the protein, and as a result the geometry (and function) of the special pair would be expected to be sensitive to protein mutation. However, these calculations were performed at the X-ray structure. We performed analogous calculations (SCF/ 3-21G) at the PCM structure and found, however, that the dimerization energy reduces to 3 kcal mol-1 for Rh. Viridis; the AM1 dimerization energy is also 3 kcal mol-1. Further, using density functional theory (SVWN with the split-valence polarized SVP basis) we calculate a dimerization energy of -12 kcal mol-1 after correction for basis-set superposition error. The difference arises as SCF-based methods do not adequately describe the pyrrole-pyrrole interactions, the SVWN result being verified in principle through MP2 calculations on model compounds.54 Hence the dimer is in fact held together by mutual attraction and, as observed, is thus expected to be robust to changes in the protein environment. For Rb. sphaeroides, our AM1-calculated dimerization energy is -1 kcal mol-1, the additional attraction arising largely from the chelation of the magnesium atom on PL by the 2a-acetyl group from PM.

Figure 5. Structure of the AM1 optimized imidazole-bacteriochlorophyll b dimer (top) compared to the PCM-optimized special pairs of Rb. sphaeroides (middle) and Rh. Viridis (bottom).

Ring Deformations in the Special Pair. The deformation of the porphyrin rings of the special pair is often explained as arising from interactions between the special pair and the protein, particularly to electrostatic forces acting between polar groups.6 To compare such influences of the surrounding protein on the dimer, a control system was investigated consisting of a bacteriochlorophyll b dimer with axial imidazole ligands in the gas phase. Its structure was optimized using AM1, constrained to maintain the pairwise ring-I interatomic distances. While unconstrained optimization would in principle be preferred, this constraint is necessary as SCF-based methods, such as AM1, describe the interring I interaction poorly: see previous section. They predict an interring separation in error by ca. +0.5 Å, and as a consequence of this, it could surreptitiously allow the porphyrin rings to relax. Our constraint prevents this, and the resulting structure for this model imidazole-bacteriochlorophyll b dimer is shown in Figure 5; there, this is compared to the corresponding structures obtained for Rb. sphaeroides and Rh. Viridis. Qualitatively, Figure 5 shows that the ring systems in the special pair are distorted, those in the model dimer are not. Clearly, the distortion does not arise from the nature of the

Bacterial Photosynthetic Reaction Center

J. Phys. Chem. B, Vol. 103, No. 23, 1999 4913

TABLE 5: The Geometrical Proximity of the Cofactors in Rh. Wiridis center-ring to center-ring distances (Å) cofactor site A

cofactor site B

PCM

1PRC

PM-ring I PL-ring I PL-ring V PM-ring V BChlM BChlL

BChlL-ring III BChlM-ring II BChlL-ring II BChlM-ring II BPhM BPhL

6.6 7.2 6.6 6.9 10.8 10.8

6.9 7.4 7.1 7.1 10.6 10.7

bacteriochlorophyll dimer internal interactions. Quantitatively, the Mg-Mg spacing is 7.87 Å, in excellent agreement with the optimized values for Rb. sphaeroides and Rh. Viridis of 7.74 and 8.31 Å, respectively. Also, the axial Mg-N bond lengths of 2.126 and 2.140 Å are similar to those calculated for the bacteria, 2.144 and 2.117 Å, as are the relative orientations of the imidazole rings. Hence the calculated gas-phase structure is very similar to those in the protein environment, the major difference being the absence of the ring deformations. This effect thus clearly must arise from interactions with the protein. Key Geometrical Proximity Effects of the Cofactors. The structural asymmetry of the photosynthetic reaction center has recently been used to explain (partly) the origin of unidirectionality of the electron transfer via the L-branch. Hasegawa and Nakatsuiji16 have approximated the transfer integrals between the L- and M-branches of Rh. Viridis by evaluating the Fock matrix element fLP,LB (where LP and LB denote the LUMOs of the bacteriochlorophylls). In order for the transfer integral fLP,LB to be large, one of the conditions that must be met is that of proximity where the atomic orbitals between two bateriochlorophylls should overlap. The molecular orbital (MO) coefficients and nodal effects also play a role in the final value of the approximated transfer integral, but here we investigate the geometrical proximity. Table 5 gives the comparisons of the key center-ring to centerring distances of the PCM-optimized geometries of the photosynthetic reaction centers with the X-ray determined distances (as used by Hasegawa and Nakatsuiji16). Three types of interactions are noted and can be described as follows: (i) Rings I of PM and PL Interacting with Rings III of BChlL and BChlM, respectiVely. The X-ray geometry center-ring to center-ring distances are 6.9 and 7.4 Å, respectively, showing that the L-branch is more compact than the M-branch. Analyzing both the MO and proximity conditions, Hasegawa and coworkers16 have concluded that the asymmetry in the contributions of rings III of BChlL and BChlM contributed to the geometrical proximity effect whereby the proximity condition overwhelms the MO coefficient condition. The PCM-optimized interring distances are somewhat shorter than the originals, but the same qualitative description of cofactor compactness remains. Hence we would expect the conclusions drawn by Hasegawa and Nakatsuiji16 to remain valid. (ii) Rings V of PM and PL Interacting with Rings II of BChlL and BChlM, RespectiVely. The X-ray center-ring to center-ring distances of rings V of the special pair with rings II of BChlL and BChlM are found to be very similar: 7.13 and 7.14 Å, respectively. Hasegawa and co-workers16 determined that the asymmetry in the contributions of rings II of BChlL and BChlM is due to the MO condition rather than the proximity condition. After optimization the PL-BChlL and PM-BChlM center-tocenter distances become 6.6 and 6.9 Å, respectively, shorter and less symmetric than the originals. This raises the possibility

TABLE 6: AMI Calculated Gas-Phase Ionization Potentials (in eV) for the Special Pair P and Its Two Halves PM and PL (All Including Axial Imidazole Ligands), as Well as for the Corresponding Imidazole-Bacteriochlorophyll b Monomer and Dimer in the Gas Phasea system

Rh. Viridis

Rb. sphaeroides

bacteriochlorophyll b

PMb PLc P

6.97 7.02 6.86

7.03 7.10 6.63

6.99 6.99 6.85

a Koopmans’ theorem60 is used throughout. b At PCM-optimized protein geometry. c At AM 1-optimized gas-phase geometry.

that there may be an increase in the proximity effect on the L-side thereby influencing the relevant electron transfer integral. (iii) The Center-to-Center Distances between BChlM and BChlL with BPhM and BPhL, RespectiVely. The center-ring to center-ring distances of BChlM-BPhM and BChlL-BPhL show little variation between the X-ray geometry and the optimized PCM geometry. Here the variation in distances are just 0.2 and 0.1 Å, respectively. Again we would expect the conclusions made by Hasegawa and co-workers16 to remain valid. Differing Protein Influence on the L and M Chains. The electrostatic and steric influences of the protein on cofactors on the L- and M-sides are quite different; these are tabulated in Table 4, as is the total QM-MM interaction energy for each cofactor. In all cases, the cofactor on the L-side is seen to fit more snugly into the protein backbone than does its M-side counterpart, the interaction energy difference being 1, 2, 16, 7, 1, 31, and 18 kcal mol-1 for P, Bchl and BPh of Rh. Viridis, and P, BChl, BPh, and Q of Rb. sphaeroides, respectively. These energies are much larger than those associated with cofactor distortion (cf. Table 1), and hence, we see that it is the proteincofactor interaction which asymmetrically distorts the structures of the cofactors. It is presumable the snug L-side proteincofactor interactions which facilitates the improved L-side intercofactor interactions mentioned above. Ionization Potentials of the Special Pair. AM1-calculated ionization potentials (IPs) are given in Table 6. Typically, these are most reliably evaluated as the difference in the energies of the neutral molecule and its radical cation, the latter being treated using spin-unrestricted Hartree-Fock (UHF) theory. UHF methodologies are, however, inappropriate for porphyrins as they incorrectly yield triplet or quartet ground states. An alternative more approximate method uses Koopmans’ theorem64 to estimate ionization energies from the orbital energies of the neutral molecule, and this is the approach which we have taken. Table 6 gives Koopmans’ IPs evaluated for imidazolebacteriochlorophyll b complex in the gas phase, evaluated both at the geometries optimized for PL and PM in Rb. sphaeroides and Rh. Viridis as well as that for the fully optimized complex. All five values are very similar, ranging from 6.97 to 7.10 eV. Given that AM1-calculated ionization potentials for magnesium compounds are usually slightly larger than experimental values,47 these values are in good general agreement with the experimental values for porphyrin related systems which usually range from 6.0 to 6.5 eV.65 Also shown in Table 6 are IPs for the PL-PM dimers as well as the model dimer from Figure 5. These indicate that dimer formation in the gas phase lowers the ionization energy, as is expected. For the model system and Rh. Viridis, this lowering is quite modest at ca. 0.14 eV, but for Rb. sphaeroides it is much larger, 0.40 eV. A very important property of the reaction center is ∆IP, the difference in the IP of PL and PM in the protein. This difference

4914 J. Phys. Chem. B, Vol. 103, No. 23, 1999 is very important in determining whether the charge in the special-pair radical cation localizes on one cofactor or delocalizes over both.25,27,33 Spectroscopic evidence66-69 places an absolute upper bound of 0.3 eV on this difference, the value being appropriate for all functioning mutants prepared to date. Although an authoritative analysis of these spectra is yet to be established, preliminary estimates25,27 suggest that |∆IP| e ca. 0.06 eV. Our calculated gas-phase values of 0.05 eV for Rh. Viridis and 0.07 eV for Rb. sphaeroides are of the expected order, as indeed is the ab initio SCF value of 0.07 eV for Rh. Viridis calculated by Sukama et al.29 In addition to indirect structural effects, the protein environment modifies ∆IP directly through asymmetry in the long-range electrostatic interactions.13,14,24,28-31 While such asymmetry is embodied in our QM/MM interaction, no effort has been made to provide a good representation. Indeed, protein models designed to account for this are very complex (see, e.g., refs 14 and 29) and, to date, are hampered by many factors such as the uncertain state of ionization of distant residues. As an example of the problems faced, Sukama et al.29 calculate that for Rh. Viridis, ∆IP changes from a gas-phase value of 0.07 eV to an in vivo value of 0.66 eV, one which is clearly an order of magnitude too large. Here we make no attempt to evaluate the effect of the protein environment on ∆IP but suggest that, in future calculations of this type, use of refined cofactor coordinates rather than the X-ray structures is essential. Conclusions The optimization of the cofactors in the bacterial photosynthetic reaction centers of Rhodopseudomonas Viridis and Rhodobacter sphaeroides using the described QM/MM method yields structures of the prosthetic groups closely preserved to the crystallographic data. Indeed, the calculated hydrogen bonding distances between GluL104 and the 9-keto group of BChlL for both protein models verifies the postulate that this amino acid residue is protonated. For PL of Rb. sphaeroides and both PL and PM of Rh. Viridis, the 2a-acetyl group points outward and is hydrogen bonded to the protein; their methyl end groups fall close to the neighboring Mg atom and pyrrole nitrogen atoms. For PM of Rb. sphaeroides, however, hydrogen bonding to the protein is restricted and these calculations predict that the 2aacetyl group rotates to coordinate its oxygen (instead of a methyl hydrogen) as a sixth ligand to the neighboring Mg, significantly increasing the interbacteriochlorophyll attraction. Deformations of the special-pair ring systems are seen not to arise from interbacteriochlorophyll interactions but rather from the much larger protein-cofactor interactions. Always, it is found that the cofactors on the L-side fit more snugly together and also into the protein than do those on the M-side. The most important result of these calculations is the production of realistic coordinates for the cofactors and the surrounding protein, ones which can safely be used in any calculation of cofactor properties. This is demonstrated most clearly in Table 1 in which extremely large differential relaxation energies are reported and in the overturning of the conclusion apparent from the work of Sakuma et al.29 that the special pair interaction is repulsive with the dimer being held together through protein compression. Acknowledgment. This work was supported by the Australian Research Council by a Postdoctoral Fellowship (M. C. Hutter and J. M. Hughes) and a Senior Research Fellowship (J. R. Reimers). We thank Dr. Charles Collyer (Department of Biochemistry) for helpful discussions concerning the interpretation of the X-ray crystallographic data.

Hutter et al. Supporting Information Available: Two ASCII data files are provided containing the optimized structures of Rh. Viridis and Rb. sphaeroides. These are in “.hin” format (rhvd.hin and rbsp.hin) for use with HYPERCHEM55 and other programs. References and Notes (1) Lendzian, F.; Lubitz, W.; Scheer, H.; Hoff, A. J.; Plato, M.; Tra¨nkle, E.; Mo¨bius, K. Chem. Phys. Lett. 1988, 148, 377. (2) Lendzian, F.; Huber, M.; Isaacson, R. A.; Endeward, B.; Plato, M.; Bo¨nigk, B.; Mo¨bius, K.; Lubitz, W.; Feher, G. Biophys. Biochim. Acta 1993, 1183, 139. (3) Mattioli, T. A.; Hoffman, A.; Robert, B.; Schrader, B.; Lutz, M. Biochemistry 1991, 30, 4648. (4) Rautter, J.; Lendzian, F.; Schulz, C.; Fetsch, A.; Kuhn, M.; Lin, X.; Williams, J. C.; Allen, J. P.; Lubitz, W. Biochemistry, 1995, 34, 8130. (5) Deisenhofer, J.; Michel, H. Science 1989, 245, 1463. (6) Ermler, U.; Fritzsch, G.; Buchanan, S. K.; Michel, H. Structure (London) 1994, 2, 925. (7) Apostolakis, J.; Muegge, I.; Ermler, U.; Fritzsch, G.; Knapp, E. W. J. Am. Chem. Soc. 1996, 118, 3743. (8) Muegge, I.; Apostolakis, J.; Ermler, U.; Fritzsch, G.; Lubitz, W.; Knapp, E. W. Biochemistry 1996, 35, 8359. (9) El-Kabbani, O.; Chang, C.-H.; Tiede, D.; Norris, J.; Schiffer, M. Biochemistry 1991, 30, 5361. (10) Boxer, S. G. Ann. ReV. Biophys. Biophys. Chem. 1990, 19, 267. (11) Treutlein H.; Schulten, K.; Bru¨nger, A. T.; Karplus, M.; Deisenhofer, J.; Michel, H. Proc. Nat. Acad. Sci. U.S.A. 1992, 89, 75. (12) Steffen, M. A.; Lao, K.; Boxer, S. G. Science 1994, 264, 810. (13) Scherer, P. O. J.; Scharnagl, C.; Fischer, S. F. Chem. Phys. 1995, 197, 333. (14) Gunner, M. R.; Nicholls, A.; Honig, B. J. Phys. Chem. 1996, 100, 4277. (15) Olson, J. M.; Miller, M.; D’Olieslager, J. Biochemistry 1995, 34, 15230. (16) Hasegawa, J.; Nakatsuji, H. J. Phys. Chem. B 1998, 102, 10420. (17) Hasegawa, J.; Ohkawa, K.; Nakatsuji, H. J. Phys. Chem. B 1998, 102, 10410. (18) Blomberg, M. R. A; Siegbahn, Per E. M.; Babcock, G. T. J. Am. Chem. Soc. 1998, 120, 18812. (19) Scherer, P. O. J.; Fischer, S. F. Chem. Phys. 1989, 131, 115. (20) Scherer, P. O. J.; Fischer, S. F. Chem. Phys. Lett. 1997, 268, 133. (21) Scherer, P. O. J.; Fischer, S. F. In The Photosynthetic Bacterial Reaction Center II: Structure, Spectroscopy, and Dynamics; Breton, J., Verme´glio, A., Eds.; Plenum Press: New York, 1992; p 193. (22) Thompson, M. A.; Zerner, M. C. J. Am. Chem. Soc. 1991, 113, 8210. (23) Thompson, M. A.; Zerner, M. C.; Fajer, J. J. Phys. Chem. 1991, 95, 5693. (24) Thompson, M. A.; Schenter, G. K. J. Phys. Chem. 1995, 99, 6374. (25) Reimers, J. R.; Hush, N. S. Chem. Phys. 1995, 197, 323. (26) Reimers, J. R.; Hush, N. S. J. Am. Chem. Soc. 1995, 117, 1302. (27) Reimers, J. R.; Hush, N. S. Chem. Phys. 1996, 208, 177. (28) Marchi, M.; Gehlen, J. N.; Chandler, D.; Newton, M. J.; J. Am. Chem. Soc. 1993, 115, 4178. (29) Sakuma, T.; Kashiwagi, H.; Takada, T.; Nakamura, H. Int. J. Quantum Chem. 1997, 61, 137. (30) Parson, W. W.; Chu, Z. T.; Warshel, A. Biochem. Biophys. Acta 1990, 1017, 251. (31) Alden, R. G.; Parson, W. W.; Chu, Z. T.; Warshel, A. J. Am. Chem. Soc. 1995, 117, 12284. (32) Ridley, J.; Zerner, M. Theor. Chim. Acta (Berlin) 1973, 32, 111. (33) Reimers, J. R.; Hutter, M. C.; Hush, N. S. Photosynth. Res. 1998, 55, 163. (34) Hutter, M. Semiempirische MO Methoden und Anwendungen fu¨ r Umgebungseffekte, Ph.D. Thesis, University of Erlangen, Germany, 1997. (35) Rauhut, G.; Alex, A.; Chandrasekhar, J.; Steinke, T.; Sauer, W.; Beck, B.; Hutter, M.; Gedeck, P.; Clark T. VAMP, Version 6.5; Erlangen, 1997. (36) Lanig, H.; Beck, B.; Clark, T. J. Mol. Model. 1995, 1, 108. (37) Lanig, H.; Beck, B.; Clark, T. Chem. Des. Automat. News 1996, 11, 19. (38) Lanig, H.; Beck, B.; Clark, T. J. Mol Model. 1996, 2, 260. (39) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227. (40) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11, 700. (41) The´ry, V.; Rinaldi, D.; Rivail, J.-L.; Maigret, B.; Ferency, G. J. Comput. Chem. 1994, 15, 269. (42) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1986, 7, 718. (43) Ranganathan, S.; Gready, J. J. Phys. Chem. B 1997, 101, 5614. (44) Bakowies, D.; Thiel, W. J. Phys. Chem. 1996, 100, 10580. (45) Gao, J.; Xia, X. Science 1992, 258, 631.

Bacterial Photosynthetic Reaction Center (46) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902. (47) Hutter, M. C.; Reimers, J. R.; Hush, N. S. J. Phys. Chem. B 1998, 102, 8080. (48) Clark T.; Gedeck, P.; Lanig, H.; Schu¨rer, G. In Molecular Modelling and Dynamics of Bioinorganic Systems; Banci, L., Comba, P., Eds.; NATO ASI Series 3, Kluwer Academic Publishers: Dordrecht, 1997; Vol. 41, p 307. (49) Rappe´, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, M. W. J. Am. Chem. Soc. 1992, 114, 10024. (50) Vasilyev, V. V.; Bliznyuk, A. A.; Voityuk, A. A. Int. J. Quantum Chem. 1992, 44, 897. (51) Dewar, M. J. S.; Thiel, W. Theor. Chim. Acta 1977, 46, 89. (52) Wang, B.; Ford, G. J. Chem. Phys. 1992, 97, 4162. (53) Bernstein, F. C.; Koetzle, T. F.; Williams, G. J. B.; Meyer, E. F.; Jr.; Brice, M. D.; Rodgers, J. R.; Kennard, O.; Shimanouchi, T.; Tasumi, M. J. Mol. Biol. 1977, 112, 535. (54) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J. Comput. Chem. 1986, 7, 230. (55) HYPERCHEM, Release 5.0 Professional Version; Hypercube, Inc.: Waterloo, Ontario, Canada, 1996. (56) Protein Data Bank at the Brookhaven National Laboratory. URL: http://www.pdb.bnl.gov. (57) Lancaster, C. R. D.; Michael, H. Photosynth. Res. 1996, 48, 65. (58) Meyer, T. E.; Donohue, T. J. In Anoxygenic Photosynthetic Bacteria, Blanckenship, R. E., Madigan, M. T., Bauer, C. E., Eds.; Kluwer Academic

J. Phys. Chem. B, Vol. 103, No. 23, 1999 4915 Publishers: Dordrecht, The Netherlands, 1995; p 725. (59) Stowell, M. H. B.; McPhillips, T. M.; Rees, D. C.; Soltis, S. M.; Abresch, E.; Feher, G. Science 1997, 276, 812. (60) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (61) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S. Karplus, M. J. Comput. Chem. 1983, 4, 187. (62) Honig, B.; Sharp, K.; Yang, A.-S. J. Phys. Chem. 1993, 97, 1101. (63) Harvey, S. C. Proteins 1989, 5, 78. (64) Koopmans, T. Physica 1933, 1, 104. (65) Depuis, P.; Roberge, R.; Sandorfy, C. Chem. Phys. Letts. 1980, 75, 434 and cited references therein. (66) Breton, J.; Nabedryk, E.; Parson, W. W. Biochemistry 1992, 31, 7503. (67) Parson, W. W.; Nabedryk, E.; Breton, J. In The Photosynthetic Bacterial Reaction Center II: Structure, Spectroscopy, and Dynamics; Breton, J., Verme´glio, A., Eds.; Plenum Press: New York, 1992; p 79. (68) Nabedryk, E.; Allen, J. P.; Taguchi, A. K. W.; Williams, J. C.; Woodbury, N. W.; Breton, J. Biochemistry 1993, 32, 13879. (69) Lancaster, C. R. D.; Michel, K.; Honig, B.; Gunner, M. R. Biophys. J. 1996, 70, 2469. (70) Clark, T.; Alex, A.; Beck, B.; Gedeck, P.; Langig, H. J. Mol. Model 1999, 5, 1.