Johnson Research Foundation, Department of Biochemistry and
Biophysics, University of Pennsylvania, Philadelphia
Pennsylvania 19104 USA
A description is given of a method to calculate the
electron transfer reorganization energy (
) in proteins using the
linear or nonlinear Poisson-Boltzmann (PB) equation. Finite difference solutions to the linear PB equation are then used to calculate
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
.
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
. 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 |
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,
G±, where
|
(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,
Go defined as the difference in
redox midpoint potentials of the donor and acceptor) and 2) the
reorganization energy,
. 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,
may be obtained by considering a pathway in which the electron is first
transferred rapidly transferred to the acceptor (t
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 
. Although calculation of
G± from
Go and
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
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
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

Go =
(Marcus and Sutin, 1985
). However,
manipulation of
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
(Marcus, 1956a
, 1963
).
To obtain an analytical expression for
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
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,
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
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
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
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
.
Nevertheless,
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
(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
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
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
in proteins.
Given the difficulty of measurement of
and the crucial role of
theory in interpreting these measurements there is still need for a
model for calculating
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
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
, 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
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
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
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
|
(2)
|
and from the diagram,
in turn is given by
|
(3)
|
where
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
G* and
Go and so should
cancel when computing
. 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
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.
The starting point is the expression for the potential distribution
(r), as a function of the position coordinate r produced by a given density distribution of charges
(r) and dipoles P(r) (Bottcher,
1973
):
|
(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,
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,
m(r). The polarization at any point is
proportional to the field and the electrostatic susceptibility
(r):
|
(5)
|
The susceptibility, defined in terms of a dielectric constant
as
(r) = (
1)/4
, can have contributions from both the electronic response (
e) and the nuclear
response (
n) so that
|
(6)
|
where
e
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:
|
(7)
|
where the shorthand notation X(
) 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
d(r) and
a(r), respectively. Substituting Eqs. 5-7
into Eq. 4, the potential distribution is given by the PB equation:
|
(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
|
(9)
|
Substituting this into Eq. 4, the potential is given by
|
(10)
|
The first and last terms cancel. Substituting for
.

d using Eq. 8 and collecting terms gives
|
(11)
|
where 
ad =
a
d
is the change in atomic charge distribution upon electron transfer, and

d*d =
d*
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:
|
(12)
|
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
|
(13)
|
where the potential in Eqs. 12 and 13 has been split into an
initial contribution,
d, and a component that changes
upon movement of charge 
x or

d*d, respectively. The fixed
contribution in Eqs. 12 and 13 is given by Eq. 8 with
x =
d.
Using Eqs. 3, 12, and 13, the reorganization energy is
|
(14)
|
Thus the reorganization energy depends only on that part of the
potential that changes upon charge transfer. It is also clear that
is invariant with respect to the direction of electron transfer;
changing the sign of dq changes the sign of 
, 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
zie/T is small throughout the
solvent. The exponential term in Eq. 7 (and thus Eq. 8) may be
linearized, leading to
|
(15)
|
Both potential terms in Eq. 14 now depend linearly on the values
of qi, so evaluating the integral yields
|
(16)
|
where 
d*d is given by Eq. 11, and

ad =
a
d is obtained
from Eq. 8 by subtraction as
|
(17)
|
It should be noted that in the linear case
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
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 |
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
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
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

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 
ad.
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
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
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
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.
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
. 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

ad when calculating
. 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 |
Tests of precision
To assess the numerical accuracy of the procedure for calculating
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
e, radii a1 and
a2, respectively, surrounded by a higher
dielectric (
p) spherical protein region, out to radii
a'1 and
a'2, respectively, separated by distance
r in a solvent of dielectric
s. The electron
is transferred between the centers of the two proteins. The analytical
approximation for
is (Marcus, 1956a
, 1963
)
|
(18)
|
Choosing a1 = 2 Å,
a2 = 3.1 Å, a'1 = a'2 = 4 Å,
e = 2,
p = 4, and
s = 80,
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 calculations.
(A) Comparison of analytical ( ) and calculated ( )
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 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
for one of the protein systems (Ru a5)-His46-myoglobin. The
variation in
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
The protein for which the most extensive calculations of
have
been performed is probably cytochrome c. Here
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
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 ( ), ions ( ),
protein ( ), and the total reorganization energy ( ); , 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
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 (
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
There is, unfortunately, no experimental data with which to
compare either the MD or FDPB calculations of
in cytochrome c.
However there are measured
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
(indicated by the horizontal error
bars on the figure, which were estimated using quoted ranges of
from the literature where available). The results suggest that the
major factors that control
- 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 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
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
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.
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
to the choice
of protein dielectric constant was examined for (Ru
a5)-His46-myoglobin, where
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
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
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 (
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
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 . Data are
for ruthenated myoglobin. Contributions are from solvent ( ), protein
( ), and the total reorganization energy ( ).
|
|
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
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 . , protein, membrane, and bulk solvent alone; ,
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
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
on distance (Liu and
Newton, 1995
). Fig. 7 shows a summary of experimental and calculated values of
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
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
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
is inversely proportional to distance,
does not hold.

View larger version (11K):
[in this window]
[in a new window]
|
FIGURE 7
Dependence of measured ( and ) and calculated
( and ) on distance. Data are for intramolecular electron
transfer in R. viridis reaction center ( and ),
ruthenated myoglobin, ruthenated cytochrome c, and ruthenated
cytochrome b5 ( and ).
|
|
 |
CONCLUSIONS |
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
by atomic detail simulations. On the other hand, the most widely used model for calculating
, 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
and the crucial role of theory in interpreting these
measurements, a method has been developed for calculating
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
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,
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
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
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
that incorporates structural detail of proteins, it
should be useful in the future for estimating
of proteins for which
measurements can't be made, for better understanding the physical
properties that control
, and in aiding in the design of synthetic
electron transfer proteins.
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.
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.