help button home button Biophys. J.
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS

This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Sharp, K. A.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Sharp, K. A.

Biophys J, March 1998, p. 1241-1250, Vol. 74, No. 3

Calculation of Electron Transfer Reorganization Energies Using the Finite Difference Poisson-Boltzmann Model

Kim A. Sharp

Johnson Research Foundation, Department of Biochemistry and Biophysics, University of Pennsylvania, Philadelphia Pennsylvania 19104 USA

    ABSTRACT
Top
Abstract
Introduction
Materials & Methods
Results & Discussion
Conclusions
References

A description is given of a method to calculate the electron transfer reorganization energy (lambda ) in proteins using the linear or nonlinear Poisson-Boltzmann (PB) equation. Finite difference solutions to the linear PB equation are then used to calculate lambda  for intramolecular electron transfer reactions in the photosynthetic reaction center from Rhodopseudomonas viridis and the ruthenated heme proteins cytochrome c, myoglobin, and cytochrome b and for intermolecular electron transfer between two cytochrome c molecules. The overall agreement with experiment is good considering both the experimental and computational difficulties in estimating lambda . The calculations show that acceptor/donor separation and position of the cofactors with respect to the protein/solvent boundary are equally important and, along with the overall polarizability of the protein, are the major determinants of lambda . In agreement with previous studies, the calculations show that the protein provides a low reorganization environment for electron transfer. Agreement with experiment is best if the protein polarizability is modeled with a low (<8) average effective dielectric constant. The effect of buried waters on the reorganization energy of the photosynthetic reaction center was examined and found to make a contribution ranging from 0.05 eV to 0.27 eV, depending on the donor/acceptor pair.

    INTRODUCTION
Top
Abstract
Introduction
Materials & Methods
Results & Discussion
Conclusions
References

Approximately one-third of all known enzymes are redox proteins; i.e., they transfer electrons as part of their function. Consequently, there is great interest in understanding how the physical properties of these proteins and their redox cofactors control the direction and rate of electron transfer. The key factors controlling the rate of electron transfer, ke, are the electronic coupling term, Hab, which describes the overlap of the donor and acceptor orbital systems and the free energy of activation, Delta G±, where
k<SUP><UP>e</UP></SUP>=<FR><NU>4&pgr;<SUP>2</SUP>H<SUP><UP>2</UP></SUP><SUB><UP>ab</UP></SUB></NU><DE>h</DE></FR> e<SUP><UP>−&Dgr;G<SUP>±</SUP>/kT</UP></SUP>, (1)
k and h are Boltzmann's and Planck's constants, respectively, and T is the temperature. For weakly coupled donor-acceptor systems the activated state for electron transfer is characterized by those nuclear configurations for which the energies of the donor (reactant) and acceptor (product) states are equal (Marcus, 1956b). Modeling this activated state is challenging as it is a nonequilibrium state the structure of which is usually unknown a priori, and which may be rarely sampled in simulations that start from a known equilibrium structure. The classical Marcus theory of electron transfer greatly simplifies this problem as it assumes that the potential surfaces for fluctuation of the nuclear positions about the equilibrium donor and acceptor states are harmonic and of the same width (Marcus and Sutin, 1985). With this linear response assumption the activation free energy depends on two parameters: 1) the driving force for electron transfer (the standard free energy of the redox reaction, Delta Go defined as the difference in redox midpoint potentials of the donor and acceptor) and 2) the reorganization energy, lambda . The latter is the energy necessary to distort the nuclear configuration from its equilibrium donor state to the acceptor state without transfer of an electron. Alternatively, lambda  may be obtained by considering a pathway in which the electron is first transferred rapidly transferred to the acceptor (t approx  1 fs) so that there is no relaxation of the nuclear configurations from their equilibrium donor state, followed by relaxation of the nuclear coordinates to the equilibrium acceptor state. The free energy change for the relaxation step is equal to -lambda . Although calculation of Delta G± from Delta Go and lambda  still involves a nonequilibrium state, the properties of this state are easier to model than those of the activated state as it involves only the equilibrium nuclear positions.

Despite the importance of lambda  in electron transfer kinetics, precise measurements of this quantity are rather rare as it is difficult to measure the reorganization energy directly. The best available method is to measure the dependence of the electron transfer rate upon Delta Go (Chang et al., 1991; Cowan et al., 1988; Gunner and Dutton, 1989; Jacobs et al., 1991; Meade et al., 1989; Therien et al., 1991; Wasielewski et al., 1985; Wuttke et al., 1992), using the fact that from Marcus theory the rate will be maximal when -Delta Go = lambda  (Marcus and Sutin, 1985). However, manipulation of Delta Go by cofactor substitution or protein mutagenesis is not possible in all systems of interest.

As reorganization occurs in response to the movement of the electron charge, it is primarily electrostatic in origin. Marcus first derived a classical electrostatics model for intermolecular electron transfer based on the Poisson equation to calculate lambda (Marcus, 1956a, 1963). To obtain an analytical expression for lambda  with this model, the donor and acceptor molecules were represented as spherical regions of radii ad and aa, respectively, separated by a distance r, and immersed in a solvent of some dielectric epsilon s. A straightforward elaboration of the model allows for the cofactors of each molecule to be surrounded by a second spherical region between it and the solvent, which in biological systems would represent the protein, with a different dielectric, epsilon p. Thus the reorganization energy has contributions from the redox cofactors themselves and from the protein(s)/solvent environment, which, for obvious reasons, became known as the inner and outer sphere reorganization energies, respectively. Although this classical electrostatics model for lambda  lacks atomic detail and it is necessarily limited by the spherical cofactor/protein assumption, it has been widely used in interpreting experimental and other theoretical studies of electron transfer as it successfully captures the general dependence of lambda  on cofactor separation, the polarizability of the medium in which the cofactors are placed, and the effect of the solvent. It is, however, only applicable to intermolecular electron transfer.

Modern simulation methods such as molecular dynamics (MD) allow one to use atomic detail models for electron transfer in protein/cofactor systems. With these methods, estimates of lambda  have now been obtained for several proteins, including cytochrome c (Muegge et al., 1997; Warshel and Churg, 1983), the photosynthetic reaction center (Marchi et al., 1993; Parson et al., 1990; Treutlein et al., 1992), and plastocyanin (Ungar et al., 1997). These simulations have greatly increased our understanding of the molecular contributions to lambda . Nevertheless, lambda  is computationally expensive to simulate by these methods, especially as it requires consideration of fluctuations from the average structure or examination of an ensemble of structures. Examination of the potential surfaces obtained from these simulations shows that, to a good approximation, they are harmonic; i.e., the linear response approximation is surprisingly accurate (King and Warshel, 1990; Marchi et al., 1993; Muegge et al., 1997; Ungar et al., 1997; Warshel and Churg, 1983). Indeed this feature can be used in the MD simulations themselves to more efficiently estimate lambda  (Muegge et al., 1997). The general applicability of linear response models for protein electrostatic phenomena is supported by other studies of the electrostatic response to proton movements (Levy et al., 1991), and the simulation of protein dipole fluctuations (Simonson and Brooks, 1996; Simonson and Perahia, 1995; Smith et al., 1993), and is no doubt one of the reasons for the continuing success of the Marcus theory of electron transfer.

Marcus's original electrostatic model for lambda  also invoked the linear response approximation, specifically in the representation of solvent reorganization by a dielectric response (Marcus, 1956a, 1963). A linear dielectric model for the solvent has also been combined successfully with self-consistent quantum mechanical calculations of lambda  for intramolecular electron transfer in small molecules (Liu and Newton, 1995). Linear response models of electrostatics, including those based on numerical solutions to the Poisson-Boltzmann (PB) equation, have also proven to be highly successful for modeling many equilibrium electrostatic properties of proteins, including redox midpoint potentials (Bashford et al., 1988; Churg and Warshel, 1986; Gunner and Honig, 1990; Langen et al., 1992). A key ingredient in this success has been the ability to rapidly and accurately solve the PB equation for essentially any arbitrary charge and dielectric distribution using finite difference (FD) and finite element (FE) methods (Holst and Saied, 1993; Nicholls and Honig, 1991; Rashin and Namboodiri, 1987; Zauhar and Morgan, 1985). This has allowed for the explicit incorporation of atomic resolution information provided by x-ray crystallography and nuclear magnetic resonance (NMR), including protein shape, positioning of charges and cofactors, and their accessibility to solvent. For many phenomena the advantages of realistically representing the structural detail of the molecule outweighs any approximations entailed by use of a linear dielectric model. Another important feature of these models is the inclusion of realistic, though implicit, solvent and solvent ion screening effects. For solvation free energy calculations, implicit solvent models often give equivalent or better accuracy compared with explicit solvent representations using MD or Monte Carlo simulations, but with one to two orders of magnitude less computation (Ewing and Lybrand, 1993; Jean-Charles et al., 1990; Simonson and Brunger, 1994; Sitkoff et al., 1994). These advances in modeling the electrostatic properties of proteins and protein/cofactor assemblies have been made possible by improvements in numerical methods and by the rapid increase in computer power. They have not yet, however, been applied to the problem of calculating lambda in proteins. Given the difficulty of measurement of lambda  and the crucial role of theory in interpreting these measurements there is still need for a model for calculating lambda  that incorporates the molecular detail provided by x-ray crystallography and NMR, yet is rapid to compute. This will facilitate the study of the physical factors controlling lambda  and the interpretation of experiment and more detailed simulations. To fill this need, in this work the finite difference Poisson-Boltzmann (FDPB) approach successfully used for equilibrium electrostatic properties is extended to the calculation of lambda , very much in the spirit of Marcus's original model, with the additional feature that the model is also applicable to intramolecular electron transfer

The outline of the paper is as follows. In the next section, Theory, a brief description is given of the nonequilibrium electrostatic theory necessary for computation of lambda  using the linear and nonlinear PB equations. In addition to showing how the contribution of solvent ion reorganization may be included, the derivation presented here is more suited to calculations with the FDPB method (and somewhat simpler) than Marcus's original derivation because the final expression for the energy does not require integrals over space of the polarization density but summations of potentials at the point atomic charges. Materials and Methods gives details of the numerical solution of the equations using the finite difference method. In Results and Discussion, the precision of the method is assessed by comparison with analytical solutions for spherical donor and acceptor molecules and by comparison with previous MD simulations on cytochrome c. The method is then applied to the photosynthetic reaction center and various ruthenated heme proteins, and comparison is made with experimental estimates of lambda  for these proteins. The effect of different physical properties such as cofactor separation and protein polarizability is also examined.

    THEORY

The standard Marcus diagram describes the energy surfaces of the donor (reactant) and acceptor (product) states as a function of the nuclear configuration coordinate by two identical parabolas, vertically offset by Delta Go, the equilibrium electron transfer free energy difference (Fig. 1). When the curvatures of the energy surfaces are the same, the energy of the activated state (where the curves cross) is given by
&Dgr;G<SUP><UP>±</UP></SUP>=(&Dgr;G<SUP><UP>o</UP></SUP>+&lgr;)<SUP>2</SUP>/4&lgr;, (2)
and from the diagram, lambda  in turn is given by
&lgr;=&Dgr;G*−&Dgr;G<SUP><UP>o</UP></SUP>, (3)
where Delta G* is the free energy change when the electron is transferred from the equilibrium donor state (d) to the acceptor state with no nuclear reorganization (d*), as indicated by the vertical transition. It should be noted that any quantum mechanical contribution to the electron transfer appears in both Delta G* and Delta Go and so should cancel when computing lambda . In other words, it is assumed, as in Marcus's original treatment, first that the cofactor contribution is small and second that the contribution of the protein and solvent to lambda  can be described via classical electrostatics using a linear dielectric response. This requires evaluation of the differences in electrostatic free energies between states d, a, and d*, where the latter is not at electrostatic equilibrium.


View larger version (14K):
[in this window]
[in a new window]
 
FIGURE 1   Marcus diagram for electron transfer. Symbols are defined in the text.

The starting point is the expression for the potential distribution phi (r), as a function of the position coordinate r produced by a given density distribution of charges rho (r) and dipoles P(r) (Bottcher, 1973):
∇ · ∇&phgr;=<UP>−</UP>4&pgr;&rgr;+4&pgr;∇ · <UP><B>P</B></UP> (4)
This equation is valid whether or not the charges and dipoles are at equilibrium. It is assumed that the potential distribution arises from some distribution of atomic charge, rho f(r), which encodes information about the structure, polarity, ionization state, and redox state of the protein/cofactor system. This distribution could represent either a time-averaged equilibrium or instantaneous nonequilibrium state of the system. The dielectric response of the cofactor, protein, and solvent gives rise to an induced dipole or polarization distribution, P(r), whereas the screening response of solvent ions gives rise to a mobile ion charge density distribution, rho m(r). The polarization at any point is proportional to the field and the electrostatic susceptibility chi (r):
<UP><B>P</B></UP>=<UP>−</UP>&khgr;∇&phgr; (5)
The susceptibility, defined in terms of a dielectric constant epsilon  as chi (r) = (epsilon  - 1)/4pi , can have contributions from both the electronic response (chi e) and the nuclear response (chi n) so that
4&pgr;&khgr;=4&pgr;(&khgr;<SUB><UP>e</UP></SUB>+&khgr;<SUB><UP>n</UP></SUB>)=(&egr;<SUB><UP>e</UP></SUB>−1)+(&egr;−&egr;<SUB><UP>e</UP></SUB>), (6)
where epsilon e approx  2 is the contribution to the dielectric constant arising from the electronic response. In the PB model the mobile charge density depends on the potential and the bulk solvent concentration of each ion, ci, of valence zi:
4&pgr;&rgr;<SUB><UP>m</UP></SUB>=4&pgr;e <LIM><OP>∑</OP><LL><UP>i</UP></LL></LIM> c<SUB><UP>i</UP></SUB>z<SUB><UP>i</UP></SUB>e<SUP><UP>−z<SUB>i</SUB>e&phgr;/kT</UP></SUP>=X(&phgr;), (7)
where the shorthand notation X(phi ) will be used to denote the functional dependence of mobile charge density on the potential.

In the initial donor and final acceptor states, the electronic, nuclear, and ionic responses are in equilibrium with their atomic charge distributions rho d(r) and rho a(r), respectively. Substituting Eqs. 5-7 into Eq. 4, the potential distribution is given by the PB equation:
∇ · &egr;∇&phgr;<SUP><UP>x</UP></SUP>=<UP>−</UP>4&pgr;&rgr;<SUP><UP>x</UP></SUP><SUB><UP>f</UP></SUB>−X(&phgr;<SUP><UP>x</UP></SUP>), (8)
where the superscript x = d or a refers to the donor or acceptor states, respectively. If the electron is moved rapidly to the acceptor, the nuclear and ionic response is fixed in the donor state, whereas the electronic polarization re-equilibrates to the field resulting from a combination of the acceptor atomic charge distribution and the donor nuclear and ionic polarization. This results in a nonequilibrium state, d*, with polarization given by
4&pgr;<UP><B>P</B></UP><SUP><UP>d</UP>*</SUP>=4&pgr;<UP><B>P</B></UP><SUP><UP>d</UP></SUP><SUB><UP>n</UP></SUB>+4&pgr;<UP><B>P</B></UP><SUP><UP>d</UP>*</SUP><SUB><UP>e</UP></SUB>=<UP>−</UP>(&egr;−&egr;<SUB><UP>e</UP></SUB>)∇&phgr;<SUP><UP>d</UP></SUP>−(&egr;<SUB><UP>e</UP></SUB>−1)∇&phgr;<SUP><UP>d</UP>*</SUP> (9)
Substituting this into Eq. 4, the potential is given by
∇ · ∇&phgr;<SUP><UP>d∗</UP></SUP>=<UP>−</UP>4&pgr;&rgr;<SUP><UP>a</UP></SUP>−X(&phgr;<SUP><UP>d</UP></SUP>)−∇ · &egr;∇&phgr;<SUP><UP>d</UP></SUP>+∇ · &egr;<SUB><UP>e</UP></SUB>∇&phgr;<SUP><UP>d</UP></SUP>−∇ · &egr;<SUB><UP>e</UP></SUB>∇&phgr;<SUP><UP>d∗</UP></SUP>+∇ · ∇&phgr;<SUP><UP>d∗</UP></SUP> (10)
The first and last terms cancel. Substituting for nabla .epsilon nabla phi d using Eq. 8 and collecting terms gives
∇ · &egr;<SUB><UP>e</UP></SUB>∇&dgr;&phgr;<SUP><UP>d∗d</UP></SUP>=<UP>−</UP>4&pgr;&dgr;&rgr;<SUP><UP>ad</UP></SUP>, (11)
where delta rho ad = rho a - rho d is the change in atomic charge distribution upon electron transfer, and delta phi d*d = phi d- phi d is the accompanying change in potential due to re-equilibration of electronic polarization only.

Assuming a representation of the charge distribution as a set of atomic point charges, qi, the electrostatic contribution to the free energy difference between initial donor and final acceptor states is given by the work performed upon moving the charge:
&Dgr;G<SUP><UP>o</UP></SUP><SUB><UP>el</UP></SUB>=<LIM><OP>∑</OP><LL><UP>i</UP></LL></LIM> <LIM><OP>∫</OP><LL><UP>q</UP><SUP><UP>d</UP></SUP><SUB><UP>i</UP></SUB></LL><UL><UP>q</UP><SUP><UP>a</UP></SUP><SUB><UP>i</UP></SUB></UL></LIM> &phgr;<SUP><UP>x</UP></SUP><SUB><UP>i</UP></SUB>(q<SUB><UP>i</UP></SUB>, q<SUB><UP>j</UP></SUB>…)dq<SUB><UP>i</UP></SUB> (12)
=<LIM><OP>∑</OP><LL><UP>i</UP></LL></LIM> <LIM><OP>∫</OP><LL><UP>q</UP><SUP><UP>d</UP></SUP><SUB><UP>i</UP></SUB></LL><UL><UP>q</UP><SUP><UP>a</UP></SUP><SUB><UP>i</UP></SUB></UL></LIM>(&phgr;<SUP><UP>d</UP></SUP><SUB><UP>i</UP></SUB>+&dgr;&phgr;<SUP><UP>x</UP></SUP><SUB><UP>i</UP></SUB>(q<SUB><UP>i</UP></SUB>, q<SUB><UP>j</UP></SUB>…))dq<SUB><UP>i</UP></SUB>,
where the potential is a function of all of the atomic charges, given by Eq. 8 for the fully relaxed state. Similarly, the electrostatic free energy change upon going from the equilibrium donor state to the intermediate state d* is
&Dgr;G<SUP>*</SUP><SUB><UP>el</UP></SUB>=<LIM><OP>∑</OP><LL><UP>i</UP></LL></LIM> <LIM><OP>∫</OP><LL><UP>q</UP><SUP><UP>d</UP></SUP><SUB><UP>i</UP></SUB></LL><UL><UP>q</UP><SUP><UP>a</UP></SUP><SUB><UP>i</UP></SUB></UL></LIM> &phgr;<SUP><UP>d∗</UP></SUP><SUB><UP>i</UP></SUB>(q<SUB><UP>i</UP></SUB>, q<SUB><UP>j</UP></SUB>…)dq<SUB><UP>i</UP></SUB> (13)
=<LIM><OP>∑</OP><LL><UP>i</UP></LL></LIM> <LIM><OP>∫</OP><LL><UP>q</UP><SUP><UP>d</UP></SUP><SUB><UP>i</UP></SUB></LL><UL><UP>q</UP><SUP><UP>a</UP></SUP><SUB><UP>i</UP></SUB></UL></LIM>(&phgr;<SUP><UP>d</UP></SUP><SUB><UP>i</UP></SUB>+&dgr;&phgr;<SUP><UP>d∗d</UP></SUP><SUB><UP>i</UP></SUB>(q<SUB><UP>i</UP></SUB>, q<SUB><UP>j</UP></SUB>…))dq<SUB><UP>i</UP></SUB>,
where the potential in Eqs. 12 and 13 has been split into an initial contribution, phi d, and a component that changes upon movement of charge delta phi x or delta phi d*d, respectively. The fixed contribution in Eqs. 12 and 13 is given by Eq. 8 with rho x = rho d.

Using Eqs. 3, 12, and 13, the reorganization energy is
&lgr;=<LIM><OP>∑</OP><LL><UP>i</UP></LL></LIM> <LIM><OP>∫</OP><LL><UP>q</UP><SUP><UP>d</UP></SUP><SUB><UP>i</UP></SUB></LL><UL><UP>q</UP><SUP><UP>a</UP></SUP><SUB><UP>i</UP></SUB></UL></LIM>(&dgr;&phgr;<SUP><UP>d∗d</UP></SUP><SUB><UP>i</UP></SUB>−&dgr;&phgr;<SUP><UP>x</UP></SUP><SUB><UP>i</UP></SUB>)dq<SUB><UP>i</UP></SUB> (14)
Thus the reorganization energy depends only on that part of the potential that changes upon charge transfer. It is also clear that lambda  is invariant with respect to the direction of electron transfer; changing the sign of dq changes the sign of delta phi , and the sign of the product is unaltered. The special case where the Poisson-Boltzmann equation can be linearized is now considered. This is valid when the protein charge density is sufficiently low and/or the mobile ion concentration is high enough that phi zie/T is small throughout the solvent. The exponential term in Eq. 7 (and thus Eq. 8) may be linearized, leading to
X(&phgr;)=8&pgr;<UP>e</UP><SUP>2</SUP>I&phgr;/kT (15)
Both potential terms in Eq. 14 now depend linearly on the values of qi, so evaluating the integral yields
&lgr;=<FR><NU>1</NU><DE>2</DE></FR> <LIM><OP>∑</OP><LL><UP>i</UP></LL></LIM>(&dgr;&phgr;<SUP><UP>d∗d</UP></SUP><SUB><UP>i</UP></SUB>−&dgr;&phgr;<SUP><UP>ad</UP></SUP><SUB><UP>i</UP></SUB>)dq<SUP><UP>ad</UP></SUP><SUB><UP>i</UP></SUB> (16)
where delta phi d*d is given by Eq. 11, and delta phi ad = phi a - phi d is obtained from Eq. 8 by subtraction as
∇ · &egr;∇&dgr;&phgr;<SUP><UP>ad</UP></SUP>=<UP>−</UP>4&pgr;&dgr;&rgr;<SUP><UP>ad</UP></SUP>−8&pgr;<UP>e</UP><SUP>2</SUP>I&dgr;&phgr;<SUP><UP>ad</UP></SUP>/kT (17)
It should be noted that in the linear case lambda  is obtained as the difference in electrostatic free energies of two equilibrium potential distributions. This point has been discussed in more detail elsewhere (Liu and Newton, 1994; Marcus, 1956a, 1963; Newton and Friedman, 1988). When the ionic atmosphere introduces a nonlinear response (Eq. 14), it is shown here that lambda  can similarly be obtained from the free energy of equilibrium potential distributions, the difference being that one must integrate over the charge transfer process.

    MATERIALS AND METHODS
Top
Abstract
Introduction
Materials & Methods
Results & Discussion
Conclusions
References

Protein coordinates

The following protein structures were obtained from the Brookhaven protein data base: the photosynthetic reaction center from Rhodopseudomonas viridis, entry 1PRC (Deisenhofer et al., 1989); sperm whale myoglobin, entry 5MBN (Takano, 1977); oxidized cytochrome b5 1CYO (Mathews et al., 1979); and oxidized tuna cytochrome c, entry 3CYT, and reduced tuna cytochrome c, entry 5CYT (Takano and Dickerson, 1981). For the homologous cytochrome c transfer, the structures of reduced and oxidized cytochrome c were positioned using INSIGHT (Molecular Simulations, San Diego, CA) with the most exposed edges of the hemes facing each other, representing the most favorable orientation for electron transfer, at center-to-center distances ranging from 18 Å (proteins in contact) to 53 Å. Ruthenated forms of myoglobin, cytochrome c, and cytochrome b5 were generated by first building the structures of the ruthenium cofactors bipyridyl, imidazole ruthenium (Ru(Bipy)Im) and pentaammine ruthenium (Ru a5) using the BUILDER module of INSIGHT. The appropriate ruthenium compound was then bonded to the most exposed nitrogen (NE2 or ND1) of the appropriate histidine and then geometry minimized while keeping the residues surrounding the histidine restrained to their x-ray crystal structure positions. For one case (myoglobin) different His48-Ru a5 positions were generated by using the BIOPOLYMER module of INSIGHT to generate the allowed rotamers of the His-46 side chain. Of the six histidine side-chain rotamers, two produced severe clashes of the side chain or the Ru a5 with the rest of the protein and were discarded. The other four rotamers were used to examine the sensitivity of lambda  to the precise cofactor position.

Calculation of potentials and reorganization energies

Electrostatic energies were calculated using the finite difference Poisson-Boltzmann (FDPB) method implemented in the software package DelPhi (Gilson et al., 1988; Nicholls et al., 1991; Sharp and Honig, 1990). Parameters used in the FDPB calculations were as follows: grid dimensions were 97 × 97 × 97, except for calculations on the photosynthetic reaction center, for which a 161 × 161 × 161 grid was used. The scale was adjusted so that the longest dimension of the protein was 80% of the grid length. The photosynthetic reaction center was placed in a large slab 40 Å thick, which was assigned a dielectric constant of 2 to simulate the electrostatic effect of the membrane. Solutions were obtained for the linear PB equation with Debye-Huckel-type boundary conditions (Gilson et al., 1988) using the multigridding method of iteration (Holst and Saied, 1993; Sharp et al., 1995), combined with dielectric smoothing and charge anti-aliasing (Bruccoleri et al., 1996) with a final convergence value of 1 × 10-4 kT/e total residual error in the potential. Radii were taken from the PARSE parameter set (Sitkoff et al., 1994). As dictated by Eq. 16, charge values were assigned to atoms in the donor and acceptor cofactors to represent the change in charge for each electron transfer being studied. The electron charge is delocalized over the cofactor atom. However, the precise distribution requires a detailed condensed phase quantum mechanical calculation. Following Gunner and Honig (Gunner and Honig, 1991), a simple formal charge scheme was adopted (Table 1), in which the charge was placed on the metal in metal-containing cofactors, on the pyrrhole nitrogens of pheophytins, or on the oxygens of quinone cofactors. To compute lambda  for a given transfer, three calculations were performed. In the first, the ionic strength was set to zero and a dielectric constant of 2 (representing electronic polarizability) was assigned to the protein and solvent to obtain delta phi d*d. In the second, dielectric constants of 4 and 80 (representing the contribution from nuclear reorientation and electronic polarizability) were assigned to the protein and solvent, respectively, and an ionic strength of 0.15 M (representing the contribution of ion reorganization) was assigned to the solvent to obtain delta phi ad. lambda  was obtained using Eq. 16 by summing over the product of the charge and the potentials. In addition, the separate contributions of protein and solvent to lambda  were estimated by repeating the calculation assigning a uniform dielectric constant of 4 and zero ionic strength to the solvent as well as the protein. The difference between these values and lambda  calculated using the full response of the solvent (dielectric constant 80, non-zero ionic strength) gives an estimate of the solvent contribution, using a donor/acceptor pair embedded in an infinitely large protein matrix as a reference. For each set of lambda  calculations, eight FDPB runs were performed with the protein mapped onto the lattice in a different position, and the average and standard deviation of the reorganization energies were computed. This translational averaging reduces the error due to the finite resolution of the lattice and provides an estimate of the numerical precision of the calculations.

                              
View this table:
[in this window]
[in a new window]
 
TABLE 1   Change in cofactor atomic charges

Treatment of crystal waters

The structure of the photosynthetic reaction center from R. viridis contains more than 200 localized waters, many of them buried (Deisenhofer et al., 1989). If there is reorientation of these buried waters upon electron transfer this could contribute significantly to lambda . The program SURFCV (Sridharan et al., 1992) was used to identify buried waters, defined as having no solvent-accessible area. To provide an upper estimate of the contribution from these buried waters, the change in electrostatic field upon electron transfer at each water site was calculated using the FDPB method, and then the waters were assumed to completely realign to this field. The change in potential at the donor and acceptor from the realigned waters was then calculated using the FDPB method, and this contribution was added to delta phi ad when calculating lambda . Surface waters seen in the crystal structure were treated like bulk solvent (i.e., they were removed from the coordinate file before calculation, so that the entire volume of the solvent exterior to the protein is assigned a dielectric of 80).

    RESULTS AND DISCUSSION
Top
Abstract
Introduction
Materials & Methods
Results & Discussion
Conclusions
References

Tests of precision

To assess the numerical accuracy of the procedure for calculating lambda  using Eq. 16, FDPB solutions were obtained for a test case for which there is an analytical approximation. The test system consists of two spherical molecules each of which has an inner cofactor sphere, dielectric epsilon e, radii a1 and a2, respectively, surrounded by a higher dielectric (epsilon p) spherical protein region, out to radii a'1 and a'2, respectively, separated by distance r in a solvent of dielectric epsilon s. The electron is transferred between the centers of the two proteins. The analytical approximation for lambda  is (Marcus, 1956a, 1963)
&lgr;=<FENCE><FR><NU>1</NU><DE>&egr;<SUP><UP>e</UP></SUP></DE></FR>−<FR><NU>1</NU><DE>&egr;<SUP><UP>p</UP></SUP></DE></FR></FENCE><FENCE><FR><NU>1</NU><DE>2a<SUB>1</SUB></DE></FR>+<FR><NU>1</NU><DE>2a<SUB>2</SUB></DE></FR></FENCE> (18)
+<FENCE><FR><NU>1</NU><DE>&egr;<SUP><UP>p</UP></SUP></DE></FR>−<FR><NU>1</NU><DE>&egr;<SUP><UP>s</UP></SUP></DE></FR></FENCE><FENCE><FR><NU>1</NU><DE>2a′<SUB>1</SUB></DE></FR>+<FR><NU>1</NU><DE>2a′<SUB>2</SUB></DE></FR></FENCE>
−<FENCE><FR><NU>1</NU><DE>&egr;<SUP><UP>e</UP></SUP></DE></FR>−<FR><NU>1</NU><DE>&egr;<SUP><UP>s</UP></SUP></DE></FR></FENCE><FR><NU>1</NU><DE>r</DE></FR>
Choosing a1 = 2 Å, a2 = 3.1 Å, a'1 = a'2 = 4 Å, epsilon e = 2, epsilon p = 4, and epsilon s = 80, lambda  was computed for various separations, as shown in Fig. 2 A. The agreement between analytical and numerical values is excellent, with a very small deviation arising from a combination of numerical error in the FDPB calculations and the approximation made in deriving Eq. 18, which neglects image charge effects.


View larger version (13K):
[in this window]
[in a new window]
 
FIGURE 2   Numerical accuracy of lambda  calculations. (A) Comparison of analytical (square ) and calculated (black-diamond ) lambda  values for electron transfer between two spherical cofactors of radii 2 Å and 3.1 Å surrounded by a spherical protein region of radius 4 Å, dielectric 4, immersed in a solvent of dielectric 80. (B) Effect of grid scale. Calculated lambda  for ruthenated myoglobin (longest dimension ~ 42 Å) as a function of grid resolution. A cubic grid of 161 × 161 × 161 points was used. Number of grids per Å is indicated on the abscissa. Error bars represent the precision estimated from eight separate calculations with different grid mappings (see Materials and Methods).

Fig. 2 B shows the effect of grid size on the calculated lambda  for one of the protein systems (Ru a5)-His46-myoglobin. The variation in lambda  when the grid spacing is reduced from 1 Å to 0.33 Å is within the numerical precision of the calculations, indicating that a sufficiently fine grid is being used for the calculations.

Comparison with molecular dynamics calculations of lambda

The protein for which the most extensive calculations of lambda  have been performed is probably cytochrome c. Here lambda  has been calculated for homologous transfer of an electron between two cytochrome c molecules using the difference in reduced and oxidized structures as revealed by x-ray crystallography (Warshel and Churg, 1983) or from MD simulations (Muegge et al., 1997). FDPB calculations were performed on the same system, and in addition, the effect of cytochrome c separation was examined (Fig. 3). The net value of lambda  varied from 16 kcal/mol to 20 kcal/mol, increasing with separation. The value at contact (r = 18 Å) was 16 kcal/mol, which is similar to the value of 9 to 15 kcal/mol estimated by Muegge et al. (1997) where the range of values in the MD simulations reflects differing restraints applied to the protein in the MD simulations used to keep the proteins at the correct distance and close to the experimental structure. The slightly smaller value seen in the MD simulations compared with the PB calculations may reflect the effect of the atomic restraints applied to the atoms, as this would have the effect of reducing the ability of the protein to reorganize somewhat.


View larger version (13K):
[in this window]
[in a new window]
 
FIGURE 3   Reorganization energy for intermolecular cytochrome c homologous transfer. Contributions are from solvent (open circle ), ions (black-triangle), protein (square ), and the total reorganization energy (black-diamond ); triangle , solvent reorganization energy in the absence of the protein. Solid line is the prediction from Marcus's model choosing a protein radius that matches the calculated value at r = 53 Å

Both MD and PB calculations give larger values than the 9 kcal/mol obtained from the calculations using the difference in x-ray structures. As pointed out by Muegge et al. (1997), this is probably because the crystal packing suppresses the structural changes associated with reorganization. Fig. 3 also shows the prediction of the Marcus model for intermolecular transfer between two spherical molecules (Eq. 18) where the radius of the protein was chosen to match the FDPB value at the largest distance. The qualitative dependence of lambda  on distance is similar as expected, as cytochrome c is reasonably spherical. Other things being equal, the reorganization energy is greater the further the electron moves because there must be greater nuclear movement to re-equilibrate with the new charge distribution. The FDPB model shows a somewhat smaller distance dependence, which arises from image charge (dielectric boundary) and shape-dependent effects included in the FDPB model.

It has been pointed out that one role of the protein is to provide a low reorganization environment to facilitate electron transfer in cases where the driving force (Delta Go) is small. Fig. 3 shows the separate contributions of protein, solvent, and solvent ions to the reorganization energy. The largest contribution is from the protein, whereas the solvent ions make the smallest contribution. The solvent makes a smaller contribution than the protein despite its greater ability to reorganize (higher dielectric) because it is further from the redox centers. This is illustrated by the top curve, calculated for the same cofactor geometry but without the protein present. Here the solvent can approach more closely, resulting in a much larger (25-32 kcal/mol) reorganization energy than with the protein present. This larger value is similar to, but somewhat smaller than, that seen in the MD simulations (37 kcal/mol) (Muegge et al., 1997). This difference may arise from two sources. First, in the PB calculations no contribution from heme flexibility has been included. This would have the effect of reducing the reorganization energy compared with the MD simulations. Second, the treatment of water is different. In the MD simulations the water was treated by a hybrid approach in which there were explicit waters in the inner region, a langevin dipole model for the water in an intermediate region, surrounded by an outer region modeled as a dielectric continuum. Muegge et al. (1997) note that "the contribution from the Langevin region to the reorganization energy is somewhat overestimated." This effect would be greatest for the calculations of the reorganization energy of the cofactor alone in water, and could account for their higher value. Overall, the agreement between the MD and PB is within what one would expect for the different methods.

Comparison with measured values of lambda

There is, unfortunately, no experimental data with which to compare either the MD or FDPB calculations of lambda  in cytochrome c. However there are measured lambda  values for intramolecular electron transfer for the photosynthetic reaction center (Farid et al., 1993; Gunner and Dutton, 1989; Lin et al., 1994; Moser et al., 1992) and various ruthenated heme proteins (Chang et al., 1991; Cowan et al., 1988; Jacobs et al., 1991; Meade et al., 1989; Therien et al., 1991; Wuttke et al., 1992). Calculations were performed on these systems, and the results compared in Fig. 4. Overall, the agreement is quite satisfactory, given the simplicity of the model, and the difficulty in measuring lambda  (indicated by the horizontal error bars on the figure, which were estimated using quoted ranges of lambda  from the literature where available). The results suggest that the major factors that control lambda - protein and solvent polarizability, cofactor separation and cofactor position with respect to the protein/solvent interface, are reasonable well represented in the FDPB model.


View larger version (13K):
[in this window]
[in a new window]
 
FIGURE 4   Comparison of the calculated value of lambda  with those measured for intramolecular electron transfer in R. viridis (Farid et al., 1993; Gunner and Dutton, 1989; Lin et al., 1994; Moser et al., 1992), ruthenated myoglobin (Cowan et al., 1988), ruthenated cytochrome c (Chang et al., 1991; Meade et al., 1989; Therien et al., 1991; Wuttke et al., 1992), and ruthenated cytochrome b5 (Jacobs et al., 1991).

Effect of ruthenium cofactor position

As the crystal structures of ruthenated heme proteins are not available, the position of the ruthenium-containing cofactor was generated by molecular modeling. This introduces an unavoidable source of uncertainty into the calculations. To assess the size of this uncertainty, calculations of lambda  for ruthenated myoglobin were performed using the four sterically permissible rotamers of the histidine to which the ruthenium cofactor is attached. The rotamers gave a variation in Ru to Fe distance of 7 Å. The results in Table 2 show that lambda  increases systematically with distance, but the uncertainty (defined as the standard deviation for the four rotamer values) is no larger than the numerical uncertainty for a single rotamer position. Thus the overall results for the four ruthenated heme proteins are judged to be rather insensitive to the precise cofactor position.

                              
View this table:
[in this window]
[in a new window]
 
TABLE 2   Sensitivity of calculated lambda  to ruthenium cofactor position

Effect of protein dielectric and bound water

In the area of protein electrostatics in general, a question of great interest is how much a protein can respond to charge movement. The treatment within the FDPB model of the protein as a dielectric material, though successful for the calculation of many properties, is somewhat of an idealization. Thus the sensitivity of lambda  to the choice of protein dielectric constant was examined for (Ru a5)-His46-myoglobin, where lambda  was calculated using values ranging from 2 to 20. Fig. 5 shows the result. When the protein dielectric constant is equal to the electronic dielectric constant (2), then there is no contribution to lambda  from the protein. Effectively there is no nuclear reorientation in the protein, and all of the effect comes from the solvent. As the protein dielectric constant is increased, the contribution from the protein increases. The contribution from the solvent simultaneously decreases because as the protein response increases, it increasingly shields the electron's field from the solvent. The combined contribution of protein and solvent increases, however, with protein dielectric constant. Overall, the best agreement between experimental and calculated values of lambda  shown in Fig. 5 occurs with a protein dielectric in the range 3-6 (data not shown). This supports the current view of the protein interior as a region of low average polarizability (defined as the ability to reorganize in response to a change in electrostatic field), although local regions may be quite polar due to the presence of preoriented dipolar groups. The proposed higher polarizability at the protein surface (epsilon  approx  20) based on calculations on pKa shift data (Antosiewicz et al., 1994) may reflect movement of charged side chains at the protein surface (Simonson and Brooks, 1996; Smith et al., 1993), but if present, these contributions do not appear to contribute significantly to lambda  for more buried electron transfer cofactors.


View larger version (14K):
[in this window]
[in a new window]
 
FIGURE 5   Effect of protein dielectric constant on lambda . Data are for ruthenated myoglobin. Contributions are from solvent (open circle ), protein (square ), and the total reorganization energy (black-square).

Protein-bound water can potentially play a large role in protein function, particularly that class that are localized due to burial in the protein interior. MD studies of electron transfer between BCl2 to BPh in the photosynthetic reaction center indicate that buried waters contribute to the energy fluctuations (Marchi et al., 1993; Treutlein et al., 1992) and may contribute to the measured difference in polarizability seen in the active and inactive sides of this molecule (Steffen et al., 1994). One hundred and thirty five such waters were identified in the R. viridis photosynthetic reaction center. Fig. 6 shows their contribution to lambda  in the FDPB model, using the full alignment approximation. (The effective definition of crystallographic waters is that they are localized positionally enough to be seen in x-ray structures. Little can be said about their orientational immobilization from the x-ray data, however). Their contribution ranges from 0.05 to 0.27 eV, depending on the donor-acceptor pair. The largest contribution from water is for the cytochrome to special pair (BCl2) transfer. This is because a substantial number of the buried waters lie between the donor and acceptor and therefore experience the largest change in field direction upon electron transfer. For all of the other donor/acceptor pairs, the waters lie more peripherally relative to the two cofactors and therefore experience a smaller change in field direction. Bearing in mind that the water contribution estimated here is an upper limit, the overall agreement between calculated and experimental values is similar with and without the water contribution.


View larger version (18K):
[in this window]
[in a new window]
 
FIGURE 6   Effect of buried water in R. viridis on calculated lambda . black-square, protein, membrane, and bulk solvent alone; square , including buried water. Donor/acceptor pairs for each transfer are indicated on the figure

Effect of electron transfer distance

In Marcus's original inner sphere/outer sphere model for intermolecular electron transfer (Eq. 18) there is a strong dependence of lambda  on the distance between the donor and acceptor. This distance dependence is borne out by studies on intermolecular electron transfer in cytochrome c here (Fig. 3) and previously. As pointed out, a two-sphere model is not applicable for intramolecular transfer although intramolecular electron transfer within a homologous series of radical anion/cations shows an inverse dependence of lambda  on distance (Liu and Newton, 1995). Fig. 7 shows a summary of experimental and calculated values of lambda  for intramolecular transfer for the four proteins studied here, plotted as a function of the center-to-center distance between donor and acceptor. Here there is no observable overall correlation of lambda  with distance. Examining the data from a single protein, the photosynthetic reaction center, where the protein size and shape is the same, there is still no observable correlation. The lack of correlation should not be interpreted as showing that lambda  is independent of distance. It does, however, show that for intramolecular electron transfer in proteins other structural factors, such as orientation of the cofactors and their relative distance from the protein/solvent interface, are equally important. Thus, for the proteins for which experimental data are available, the Marcus-type behavior, where lambda  is inversely proportional to distance, does not hold.


View larger version (11K):
[in this window]
[in a new window]
 
FIGURE 7   Dependence of measured (black-square and bullet ) and calculated (square  and open circle ) lambda  on distance. Data are for intramolecular electron transfer in R. viridis reaction center (open circle  and bullet ), ruthenated myoglobin, ruthenated cytochrome c, and ruthenated cytochrome b5 (square  and black-square).

    CONCLUSIONS
Top
Abstract
Introduction
Materials & Methods
Results & Discussion
Conclusions
References

The reorganization energy is a key parameter controlling the electron transfer rate in proteins. It is, however, difficult to measure the reorganization energy directly. It is also computationally expensive to simulate lambda  by atomic detail simulations. On the other hand, the most widely used model for calculating lambda , Marcus's electrostatic model, has been unnecessarily limited by the assumption that the cofactor/protein is spherical and is applicable only to intermolecular electron transfer. Given the difficulty of measurement of lambda  and the crucial role of theory in interpreting these measurements, a method has been developed for calculating lambda  that is applicable to both inter- and intramolecular transfer, which incorporates the molecular detail provided by x-ray crystallography and NMR and yet is easier to compute. The model is based on the FDPB approach successfully used for modeling equilibrium electrostatic properties of proteins.

The results of calculations on 12 electron transfer reactions in four different proteins show that the major factors controlling lambda  are captured by the FDPB model: the distance between donor and acceptor, the position of the cofactors with respect to the protein-solvent boundary, and the polarizability of the protein and solvent. For intermolecular electron transfer, lambda  is primarily dependent on the distance between donor and acceptor in a very similar way to that predicted by the Marcus model. For intramolecular electron transfer, however, there is no simple dependence on distance and the cofactors' positions with respect to the protein-solvent surface are at least as important. The major contribution to lambda  comes from the protein polarizability, with a smaller contribution from solvent and, at least in the proteins examined here, a very small contribution from solvent ions. The protein contribution is large because it surrounds the cofactors and therefore feels the greatest change in field upon electron transfer. Nevertheless, the results confirm the idea of proteins as providing a low reorganization energy environment for electron transfer as the reorganization energy for an equivalent donor/acceptor pair in solvent in the absence of the protein would be much higher. The effect of buried water in the photosynthetic reaction center is generally quite small (<0.2 eV). The results show that a significant contribution comes only when a substantial number of waters lie between the donor and acceptor. Based on these calculations, and the dependence of lambda  on the key protein structural parameters, it is possible to put a reasonable upper limit of 1.6 eV on the reorganization energy of proteins.

As the FDPB model provides a computationally tractable method for calculating lambda  that incorporates structural detail of proteins, it should be useful in the future for estimating lambda  of proteins for which measurements can't be made, for better understanding the physical properties that control lambda , and in aiding in the design of synthetic electron transfer proteins.

    ACKNOWLEDGMENTS

I thank Les Dutton, Chris Moser, and Marilyn Gunner for many helpful discussions.

This work was supported by National Science Foundation grant MCB95-06900 and the E. R. Johnson Research Foundation.

    FOOTNOTES

Received for publication 11 August 1997 and in final form 2 December 1997.

Address reprint requests to Dr. Kim A. Sharp, Department of Biochemistry/Biophysics, University of Pennsylvania, 37th & Hamilton Walk, Philadelphia, PA 19104-6059. Tel.: 215-573-3506; Fax: 215-898-4217; E-mail: sharp{at}crystal.med.upenn.edu.

    REFERENCES
Top
Abstract
Introduction
Materials & Methods
Results & Discussion
Conclusions
References