Theoretical determination of the standard reduction potential of plastocyanin in vitro

Article abstract of DOI:10.1021/jp0307452

The determination of the reduction potential of the protein plastocyanin was discussed using quantum chemical QM/MM calculations. Protein was divided into smaller parts using ONIOM method and it was based on its relative contribution to reduction potentia

Full text of DOI:10.1021/jp0307452

J. Phys. Chem. B 2004, 108, 8007-8016
8007
Theoretical Determination of the Standard Reduction Potential of Plastocyanin in Vitro
Sambhu N. Datta,* Jawahar Sudhamsu, and Anshu Pandey
Department of Chemistry, Indian Institute of Technology, Bombay, Powai, Mumbai 400076, India
ReceiVed: June 11, 2003; In Final Form: March 3, 2004
Quantum chemical QM/MM calculations have been performed on the copper-containing blue protein
plastocyanin that is involved in the photosynthetic electron transfer. Crystallographic coordinates of the non-
hydrogen atoms in the oxidized and reduced forms of poplar plastocyanin were obtained from the Protein
Data Bank. The present work required calculations on the oxidized form that has a molecular structure
independent of pH, and the reduced form with different structures at pH values 3.8 and 7. At pH 7, both the
oxidized and reduced forms of the protein have distorted tetrahedral geometry for the copper-containing active
site, the Cu atom being coordinated to two histidine, one cysteine, and one methionine residues. At pH 3.8,
the active site of the reduced species is trigonally coordinated to one histidine, one cysteine, and one methionine
residues. To optimize the geometry of the system while retaining the constraints of the protein backbone, the
ONIOM methodology was adopted so as to treat the active site by the DFT-B3LYP method using the 6-31G
basis set, whereas the geometries of the nearby residues as well as six neighboring water molecules were
optimized by the MM/UFF method. Then atomic charges for the atoms of the protein (apart from those in the
active site) were determined from DFT calculations separately for each amino acid residue using the STO-
3G basis set. The atomic charges of the water molecules were computed by the DFT/6-31G(d) method. Finally,
the electronic energies were recalculated by the ONIOM technique where the DFT-optimized active site was
again treated at the 6-31G level of theory, whereas the rest of the protein, along with the solvent molecules
near the active site, were treated by the Amber force field method using the calculated DFT charge on each
atom. This treatment effectively allowed us to retain the steric constraints offered by the protein backbone
during the optimization process, as well as the effects arising from the interaction of the protein dipoles and
the bare charges on the protein with the atomic charges in the complex, thereby accommodating all electronic
interactions such as charge-charge, charge-dipole, dipole-dipole, etc. The thermal energies of various oxidized
and reduced forms were computed for a slightly simplified model of the active site (the part optimized in
ONIOM) by using the DFT-B3LYP methodology. The effective radius of the globular protein plastocyanin
was determined from the crystallographic data. The stability of each species arising from its interaction with
medium was determined by explicitly calculating the Born charge-dielectric (water) interaction energy and,
for the solvated proton, the Debye-Hu¨ckel energy of ion-ionic atmosphere interaction. The dielectric constant
of water and plastocyanin were taken as 78.5 and 8.0, respectively. The interaction with the medium and the
entropy changes are found to play a critical role in determining the reduction potential. The process of reduction
of plastocyanin in an aqueous medium involves a very large reorganization of water molecules, and a large
entropy change that cannot be computed readily. Hence, the entropy of reduction of plastocyanin was taken
from experimental data that are available for pH 7. The free energy change was calculated for the reduction
of plastocyanin in water and proton in water. From these values, the standard reduction potential was determined
at pH 7. The calculated potential (376 ( 38 mV) is in excellent agreement with the observed one (379 mV)
for the radius of the globular protein 22.37 ( 0.16 Å. A similar calculation leads us to predict the entropy of
reduction of plastocyanin at pH 3.8.
1. Introduction
Plastocyanin from poplar leaves is a protein with a single
polypeptide chain, has a molecular weight around 10 500 Da,
Plastocyanin (PC) is a 97-99 residue redox active protein1
that stays in the aqueous phase of chloroplast, near the thylakoid
membrane. It is the redox partner of cytochrome f and P700 in
the Z scheme of photosynthesis. It accepts an electron from
cytochrome f that is embedded in the condensed phase of the
thylakoid membrane and then migrates through the stroma to
reach Photosystem I, that is, P700. Spectroscopic studies on
the soluble type 1 blue copper protein have shown an intense
absorption band near 600 nm and a midpoint reduction potential
of 0.379 V at 25 °C at pH 7.2
and resembles a slightly flattened spheroid. The crystal structure
of poplar plastocyanin in its oxidized form and that of the
reduced forms at various pH have been reported.3 Recent
structure analyses of the oxidized form of plastocyanin at
different pH values showed that the oxidized Cu site was, within
the limits of precision, invariant in the range 4.2 e pH e 6.0.4
The most significant differences among the structures of the
reduced protein at different pH values, or between them and
the structure of the oxidized protein, are concentrated at the Cu
site. In both the oxidized and reduced forms at high pH, the Cu
atom is coordinated by the Nδ (imidazole) atoms of His37 and
* Corresponding author. E-mail: sndatta@chem.iitb.ac.in.
10.1021/jp0307452 CCC: $27.50 © 2004 American Chemical Society
Published on Web 05/18/2004

8008 J. Phys. Chem. B, Vol. 108, No. 23, 2004
Datta et al.
Figure 1. Active site in plastocyanin, with the binding residues
modeled as amino acids. The stereoscopic projection is shown here.
His87, the Sγ (thiolate) atom of Cys84, and the Sδ (thioether)
atom of Met92 (Figure 1). The distorted tetrahedral geometry
and the unusually long Cu-S (Met92) bond are retained in both
these forms. The only effects of the reduction at a high pH are
a lengthening of the two Cu-N (His) bonds by about 0.1 Å,
and small changes in two bond angles involving the Cu-S (Cys)
bond. As the pH is lowered, the Cu atom and the four Cu-
binding side chains appear to undergo small but concerted
movements in relation to the rest of the molecule. At a low pH
(pH 3.8), the reduced form has the Cu atom trigonally
coordinated by Nδ (His37), Sγ (Cys84), and Sδ (Met92). The
fourth Cu-ligand bond is broken, the Cu atom making only a
van der Waals contact with the imidazole ring of His87. The
imidazole ring is rotated by 180° about the Câ-Cγ bond with
respect to its orientation in the high-pH form and is protonated.5
It is known that plastocyanin has a pH-dependent reduction
potential below pH 5.5.6 Plastocyanins from different sources
show somewhat varying redox potentials. For example, whereas
poplar plastocyanin has a redox potential of 441 mV at pH 3.6,
plastocyanins from cucumber, spinach, and parsley show slightly
different values of 458, 473, and 497 mV, respectively, at the
same pH. Similar trends are noticed for other pH values.6c
Reference 6c shows that at pH 7, plastocyanin from Populus
nigra variant italica has a midpoint potential of 0.379 V. The
rate of reaction between the reduced form of plastocyanin and
inorganic oxidants depends on the pH value, and the existence
of a protonated redox-inactive form of the protein was proposed
earlier.7 It is now contended that in aqueous solution and below
pH 5.5, the redox-active high-pH form and the redox-inactive
low-pH form must be in a dynamic equilibrium. Also, the fact
that the structural changes related to the oxidation state or pH
must be highly localized is supported by independent evidence
from 13C and 1H nuclear magnetic resonance measurements on
the oxidized and reduced forms of plastocyanin in solution.8
Plastocyanins from different sources have been electrochemi-
cally characterized, and their reduction potentials at different
pH are known.6
Figure 2. Charged surface residues are shown in stick model along
with the active site. The rest of the protein is in the wire model.
the charged residues in the stick model and the rest of the protein
in the wire model. Poplar plastocyanin has seven aspartic acids,
eight glutamic acids, and six lysine residues but no arginine
residue. The two histidine residues are bound to the metal and
hence they are not protonated at the physiological pH.
As a part of our ongoing investigations on the redox processes
occurring in the Z-scheme,9-13 we have theoretically determined
the reduction potential of plastocyanin by using the quantum
mechanics/molecular mechanics (QM/MM) method. The mo-
lecular structure selected for our quantum chemical calculation
is discussed in section 2. A short discussion on the method of
calculation is given in section 3. The computed results are
presented in section 4. The energies and the molecular orbitals
computed in this work will be used to determine the rates of
electron transfer from cytochrome f to plastocyanin and finally
from plastocyanin to P700. Our conclusions are summarized in
section 5, where we also discuss the limitations of the present
work in detail.
2. Molecular Geometry
Crystallographic Geometry. The basic aim of this work is
to calculate the reduction potential of plastocyanin in vitro.
Whereas the crystallographic structure represents the structure
of the biomolecule in a condensed phase, it is also known to be
very unreliable for proteins in solution, and thus it becomes
necessary to optimize the geometry of the protein model used
for computations. A significant problem involved is that of
retention of the protein backbone during the optimization process
for the geometry of the active site and its immediate environs.
To effectively retain the geometry constraints, we used an
ONIOM methodology to simulate the nearby residues of the
protein at a low level of theory, namely, the molecular
mechanics/universal force field (MM/UFF) method, while
modeling the active site at a high level of theory, namely, density
functional treatment (DFT). The starting points of the optimiza-
tions were the protein crystallographic coordinates at different
pH values. Crystallographic coordinates were obtained from the
Protein Data Bank at Brookhaven National Laboratory (PDB
Code: 1PLC of the plastocyanin in oxidized form at pH 6.0,
5PCY of plastocyanin in reduced form at pH 7.0, and 6PCY of
plastocyanin in reduced form at pH 3.8 from poplar leaves, P.
nigra variant italica).
One important feature of the protein is that in solution, it is
highly charged in the pH range 2-10. In general, the amino
group in the side chain of lysine and arginine residues in a
protein is protonated at pH e 10 if it is exposed to the solvent.
For histidine residues the δ nitrogen atom is protonated at pH
e 6. Similarly, the carboxylate groups in the side chains of
solvent-exposed glutamic acid and aspartic acid residues are
deprotonated at pH g 2.0. Figure 2 shows the active site and

Standard Reduction Potential of Plastocyanin
J. Phys. Chem. B, Vol. 108, No. 23, 2004 8009
The copper-containing active site is close to the surface of
the molecule and is coordinated to a surface exposed histidinyl
(His-87) imidazole moiety. The three other coordinating atoms
are a cysteinyl (Cys-84) thiolate, a methioninyl (Met-92) sulfur
and another histidinyl (His-37) imidazole, as shown in Figure
1. The geometry is a distorted tetrahedron. This distortion is
thought to be responsible for the relatively high midpoint
potential of the protein. More recently, the abnormally long
copper-thioether axial bond has also been cited as a key factor
responsible for the relatively high midpoint potential of the
protein.14 The free energy of reduction must be predominantly
determined by the electronic structure of the Cu complex, which
is the redox center. The same free energy, nevertheless, is a
small quantity and so its magnitude is greatly influenced by
the surrounding protein chain. The interaction between the
complex and the enveloping proteinous part is mainly electro-
static in nature. Hence the calculation of the electronic structure
of the complexes, though exhaustively accounting for the
interaction with the medium, is expected to give rise to a
reasonably good estimate of the reduction potential.
Model Complexes. The present work required calculations
on the oxidized form that has a structure independent of pH,
and the reduced form with different structures at pH values 3.8
and 7.
The structure of the oxidized form was reported to be
essentially invariant in the pH range 4.2 e 6.6 In this work we
assume that the structure of the oxidized form also holds for
the pH range 3.8-7, in which we are computing the reduction
potentials.
The contribution of the protein backbone to the evolution of
the geometry of the active site was simulated at first by including
the nearby residues in the optimization process. Thus, the
residues proline 36, aspargine 38, tyrosine 83, serine 85, proline
86, glutamine 88, glycine 91, and valine 93, as well as a number
of solvent molecules (six) were included in the model and treated
at a low level of theory, namely, the molecular mechanics/
universal force field (MM/UFF) method, during the optimiza-
tion. The active site containing the Cu atom, the histidine
imidazole rings, and the cysteine -S-CH2 and methionine
CH3-S-CH2 groups (modeled in ONIOM as -S-CH3 and
S(CH3)2, respectively) as shown in Figure 3 were studied by
the DFT/B3LYP methodology. The neighbors and the active
site for the first set of ONIOM calculations are shown in Figure
4 for the oxidized species. The DFT atomic charges on the
protein were computed separately for each amino acid residue
using the STO-3G basis set, with the individual amino acids
being modeled as amino aldehydes, as suggested by the
reviewers. The N terminal of the amide linkages were clipped
to amines, whereas the carboxylate terminals were treated as
aldehydes. Similarly, atomic charges of the water molecules
were computed by the DFT/6-31G(d) method. The electronic
energies were determined by treating the DFT-optimized region
at a 6-31G level of theory, but the rest of the protein, along
with the 44 solvent molecules near the active site (which had
been previously optimized by Amber force field), were now
treated as point charges. This procedure not only keeps up the
steric constraints offered by the protein backbone but also
accommodates all electronic interactions such as charge-charge,
charge-dipole, dipole-dipole, etc. The final model complex
for the oxidized form is shown in Figure 5.
used to compute electronic coupling interactions in long-distance
electron transfer processes such as those involved in photosyn-
thesis. The calculation of the redox potentials of biomolecules
that are generally large species is a notoriously difficult job. In
the theoretical studies of large systems, smaller models of the
real systems have been traditionally employed. In the case of a
system as large as a protein, the underlying problem is to retain
the effects of the large protein chain on the active site while
simulating the active site at a high level of theory. The ONIOM
method offers an attractive technique to effectively study the
characteristics of a large biological system, by dividing the
system into components or layers, which are studied at differing
levels of theory. The ONIOM33 method expresses the energy
of the system as
n
E(ONIOMn) ) E[Level(i),Model(n+1-i)] -
∑
i)1
n
E[Level(j),Model(n+2-j)] (1)
∑
j)1
During geometry optimization of the active site, we have
simulated the system by a two-layer ONIOM, in which the
active site is studied using density functional theory, whereas
the lower layer consisting of adjoining residues is treated by
using molecular mechanics. Density functional theory (DFT)
has revolutionized studies on transition metal complexes and
proteins (especially blue copper proteins), providing much more
accurate results at approximately the same cost as Hartree-
Fock theory. The DFT-B3LYP/6-31G method is adopted for
the higher layer, as the species of interest involves a transition
metal ion. In the course of optimizing the geometry of the active
site with a few nearby residues, the lower layer has been treated
by the UFF method.34 The UFF method uses a charge equilibra-
tion scheme (QEq),35 which enables charges to be generated
dynamically in response to the environment and has been shown
to reproduce atomic charges accurately for a variety of systems.
Gaussian 98 on Windows (G98W)36 was used in all these
calculations. The effect of the interaction of protein dipoles and
bare charges on the protein with the active site has been directly
included into electronic energies by use of the two-layer ONIOM
scheme (Figure 5). Unfortunately, the GAUSSIAN-98 software
could not handle the whole of the protein. Therefore, the
remainder of the protein was optimized by an Amber force
field,37 with DFT charges being used on the atoms. The Amber
force field parameters generate accurate properties for large
biomolecules such as proteins, nucleic acids, etc. Additionally,
the Amber force field employs a dielectric constant for scaling
down long-range coulomb interactions, which makes it ideal
for optimizing the protein bulk that does not have any hetroatom.
The electrostatic effects of the protein on the active site were
now determined. The active site has once again been treated at
a 6-31G level of theory. The protein is assumed to have an
effectiVe dielectric constant (ꢀp), this being a good approximation
as the protein interior is fairly hydrophobic.
Zhan et al.38a have determined the absolute hydration free
energy of a proton in water. In our calculations, we make use
of the value of hydration energy obtained by Zhan et al. to
calculate the potential of the normal hydrogen electrode. The
other species at the electrode, namely, molecular hydrogen, was
optimized by using the DFT-B3LYP procedure, using the
6-31G(d) basis set. All molecular characteristics such as thermal
energy and entropy were also calculated at this level of theory.
Zhan et al. have also determined the hydration free energy of
3. Method of Calculation
Methodology. For more than two decades, ab initio,15-21
semiempirical,22-25 and even empirical26-32 methods have been

8010 J. Phys. Chem. B, Vol. 108, No. 23, 2004
Datta et al.
Figure 3. Systems chosen for thermal energy determination for the oxidized species (top) and the reduced species at pH values 7.0 (middle) and
3.8 (bottom).
the electron through a first principles approach.38b They have
assigned a value of -35.5 kcal/mol to the free energy of electron
hydration. This neglects the additional stabilization gained by
an electron in a solution of an electrolyte due to formation of
an ion cloud around the electron in accordance with Debye-
Hu¨ckel theory. The hydration energy of the electron has been
explicitly added, although the same energy has no direct bearing
on the determination of reduction potential.
Systems representing the active sites and used to estimate
the thermal energy and entropy for the species are shown in
Figure 3. These calculations have also been carried out at the
DFT/6-31G level. The thermal energy correction to the elec-
tronic energy incorporates the zero point energies for the system.
Medium Effects. The additional stability attained by a
solvated species can be considered most efficiently by consider-
ing a dielectric continuum around it, and evaluating its inter-
action by using a self-consistent reaction field method.39
However, the size of the protein is prohibitively large and does
not permit such a calculation. Thus, the solvent medium, in this
case water, has been treated as a dielectric continuum, and the
energy of interaction of the protein with the solvent has been
considered as the Born energy of ion-dielectric interaction and
the Onsager energy of dipole-dielectric interaction.
Olsson and Ryde40 performed extensive calculations on a
model of the active site, Cu(Im)2(SCH3)(S(CH3)2), where the
histidines were replaced by imidazole moieties, the cysteine was

Standard Reduction Potential of Plastocyanin
J. Phys. Chem. B, Vol. 108, No. 23, 2004 8011
In contrast, we are interested in computing the absolute value
of the reduction potential of the whole protein in vitro and not
the relative value for the active site dissolved in water. As
discussed earlier, a direct SCRF treatment of the protein is not
possible due to the prohibitively large size of the species. Also,
because the active site is buried inside the protein and exposed
to the solvent only on one side, and in turn the protein remains
in solution, a direct SCRF or a PCM calculation would not give
a true representation of the system. Any meaningful calculation
would require an extensive change in the source code. Therefore,
we have explicitly considered the Born energy as a first-order,
perturbative effect.
Torres et al.42a have performed density functional studies on
the redox potentials of Fe4S4 clusters. DFT optimized geometries
were used to undertake redox potential calculations on the five
known oxidation states of [Fe4S4(SCH3)4]n-, where n lies
between 0 and 4. The authors assumed the environment to be
a continuum dielectric during the course of optimizations. The
authors then used a standard set of radii and atomic charges to
simulate the effect of protein dipoles on the active site. A triple-ú
basis set was used for 4s and 3d valence orbitals. A single 4p
polarization function and 3s and 3p orbitals were used for iron,
whereas 3s and 3p functions used on S atoms were augmented
with a 3d polarization function.
Figure 4. Complex chosen for the first set of ONIOM optimizations
of the active site. The DFT layer is in dark, and the UFF layer is lightly
shaded.
Olsson et al.42b have recently investigated a frozen density
functional free energy simulation of redox potentials of plas-
tocyanin and rusticyanin. These authors have divided the protein
into regions, which have been treated at differing levels of
theory. Thus, the active site has been treated by the DFT
methodology. The 63321/5211*/41+ basis set had been used
for copper, the 73111/6111/1* basis set for sulfur, the 5211/
411/1 basis set for oxygen, nitrogen, and carbon, and the 41/1
basis set for hydrogen. Other protein residues have been treated
by the frozen DFT methodology, which assumes the electronic
density to be frozen in space while the nuclei are free to move.
For the fragments treated by FDFT, the smaller 6321/521 basis
set was adopted for sulfur, and the 521/41 basis for oxygen,
nitrogen, and carbon. The 41 basis set was used for hydrogen.
The methodology and the basis sets used completely differ from
those used in the present work. Two significant assumptions
made by Olsson et al. are as follows: (i) The Cu-S bond length
for the copper-cysteine thiolate linkage in the crystallographic
geometry was retained in their calculations. (ii) The dielectric
constant of the protein was taken to be 60 for the calculation of
the charge-charge interaction energy between the ionized
residues and the active site. The physical interpretation for such
a high value of dielectric constant draws from the reorientation
of charged residues of the protein.42c,d However, the protein
plastocyanin has a hydrophobic interior, and hence a high degree
of unfolding of the protein or reorientation of charged residues
may be questionable. Indeed, plastocyanin is known to be a
globular protein. Thus, the main emphasis of the work in ref
42b is to reproduce the differences in the reduction potentials
of the two blue copper proteins plastocyanin and rusticyanin.
Olsson et al. thus infer that the FDFT formalism is applicable
to the determination of electronic properties of large systems.
Figure 5. System chosen for the electronic energy calculations. The
DFT layer is dark, and the Amber layer is lightly shaded.
represented as a methyl thiolate, and the methionine residue
was modeled as dimethyl thioether. These authors suggested
that the surrounding protein dipoles have an effect on the
reduction potential. Solvation energy was calculated by the PCM
method,41 assuming a constant dielectric continuum of water
around the complex. Olsson and Ryde calculated the relative
reduction potentials for a number of Cu-containing proteins as
differences of single point energies in the solvent (water),
performed on the vacuum equilibrium geometries of the reduced
and oxidized complexes. In their calculations, the reduction
potential of plastocyanin was assumed at the experimental value,
and relative to this value, the values of reduction potentials of
other blue-copper proteins were determined.
The emphasis of the present work is to determine the value
of the reduction potential of plastocyanin relative to the standard
hydrogen electrode, as well as to determine its entropy of
reduction at low pH values from its electrochemical properties.
In our work, we have not imposed restrictions on the optimiza-
tion of the copper-cysteine thiolate linkage. Molecular dynam-
ics simulations have been undertaken to establish the appropriate
dielectric constant for the protein in an aqueous solution. The

8012 J. Phys. Chem. B, Vol. 108, No. 23, 2004
Datta et al.
results of our molecular dynamical studies are discussed in detail
in the appendix.
Ion-Dielectric Interaction. The Born energy of the protein
in water is written as
Free Energy. For the smaller species involved in the reaction
at the normal hydrogen electrode, the total internal energy is
written as
Etotal ) ESCRF + EDH + EThermal
(5)
2
zP2 eo
1
ꢀw
where EDH is the Debye-Hu¨ckel energy of interaction of the
EBorn(plastocyanin in water) ) -
1 -
(2)
(
)
2rg
ion with its ionic atmosphere,44
zi2 eo κ
2
In the above, eo is the electronic charge, ꢀw is the dielectric
constant of water, zpeo is the net charge of plastocyanin (zp is
-8 for the oxidized form and for the reduced form at pH 3.8,
and -9 for the reduced form at pH 7), and rg is the effective
radius for the protein plastocyanin, respectively.
EDH ) -
(6)
2ꢀw
The Debye-Hu¨ckel reciprocal length κ is given by
The effective radius of plastocyanin was calculated from the
crystallographic data. The dimensions of a unit cell are 29.6 Å
× 46.86 Å × 57.6 Å (PDB code 1PLC). The volume of the
molecule was found from the volume of the unit cell by
considering the packing fraction π/6, that is,
2
1/2
4πeo
Nci
κ )
zi2
(7)
∑
(
)
[
]
ꢀwkBT
1000
i
In the above, kB is the Boltzmann constant, T is the temperature
of the system in Kelvin, N is the Avogadro number, ci is the
concentration of the ith species in molar unit, and zieo is the
charge of the ith ion. The energy EDH is zero for water but finite
for the aquated proton. No medium interaction energy is
considered for the gaseous H2.
4
3
π
πrg3 ) abc
6
(3)
where a, b, and c are the side lengths of the unit cell. The value
of rg turns out to be 21.53 Å or 40.69 au (1 au ) 0.529 177 Å).
The solvent water molecules placed around the active site are
additionally found to contribute 7765 au to the volume of the
protein. Thus, the effective radius of the hydrated protein
becomes 21.73 Å or 41.06 au. The normal practice, as in
GAUSSIAN, is to add half the radius of the solvent molecule
(water) to the molecular radius. Because there are highly polar
regions in the exterior of the protein and these sites are
The total energy of the protein is written as
Etotal ) EONIOM2 + EBorn + EOnsager + EThermal
(8)
Here, EDH was not considered because of the very low
concentration of all ionic species.9c
The standard enthalpy H° and the standard Gibbs free energy
G° are given by
2
1
considerably hydrated, we decided to add /3rw ( /6rw, where
rw ) 0.96 Å. Thus the effective radius rg is found to be 22.37
( 0.16 Å or 42.27 ( 0.30 au.
Ho ) Etotal + PV
The inclusion of an adequate number (44) of water molecules
about the active site of plastocyanin in the computational model
effectively simulates the solvation of the imidazole moiety at
pH 3.8. The dipole moments (µ) of the oxidized species, reduced
species at pH 7, and reduced species at pH 3.8 are 73.00, 75.96,
and 65.20 D, respectively. These values are obtained by
summing over all the bond dipoles. The Onsager energy of
interaction between the protein dipole and the solvent, given
by
and
G° ) H° - TS°
(9)
Reductive Processes. We calculate the free energy for the
overall process
1
H+(aq) + e- (aq) f H2(g)
(10)
2
Thus, for the smaller species involved in the reaction at the
normal hydrogen electrode, the standard Gibbs free energy
change is given by
µ2
ꢀw - 1
2ꢀw + 1
EOnsager ) -
(4)
3
[
]
rg
1
2
∆G° H+(aq)/ H2(g) )
can be estimated from the computed data. The Onsager term
introduces a further correction of 12.0 mV to the reduction
potential at pH 7, and a correction of -29.5 mV at pH 3.8. The
estimate for the reduction potential at pH 7 after inclusion of
the Onsager energy is thus 376 mV.
[
]
∆ESCRF/DFT + ∆EDH + ∆Ethermal + kBT∆n - T∆S° (11)
per molecule, where ∆n is the change in the number of species
in the gas phase. This is equivalent to
Dielectric Constant of the Interior of the Protein. Molecular
dynamics studies on metal proteins such as myoglobin and
cytochrome f have shown that proteins can be better considered
as concentric shells of different effective dielectric constants,
rather than one uniform dielectric constant.43 Assuming that the
protein is a sphere of two layers of concentric dielectric media,
one for the highly polar surface (which makes the protein water
soluble) and one for the hydrophobic interior, an expression
was derived for the effective dielectric constant of the interior
of the protein. The explicit derivation is given in Appendix I.
The calculated value of ꢀp, 8, was used in the calculation of the
energy of electrostatic interaction between the distant residues
and the protein active site.
1
2
∆G° H+(aq)/ H2(g) )
[
]
1
2
-∆GHydration(proton) + ∆EDH + ∆G H+(g)/ H2(g) (12)
[
]
The reduction of PC+(aq) can be considered to occur as follows.
Above pH 5.5, the reaction is just the addition of an electron to
the Cu atom,
PC+(aq) + e-(aq) f PC(aq)
whereas below pH 5.5, the reaction that takes place is the
(13)

Standard Reduction Potential of Plastocyanin
J. Phys. Chem. B, Vol. 108, No. 23, 2004 8013
TABLE 1: Calculation of the Gibbs Free Energy for the Hydrogen Smaller Speciesa
electronic
energy (au)
thermal
Debye- Hu¨ckel
correction (au)
total
entropy
free
energy (au)
species
energyb (au)
enthalpyc (au)
(cal mol-1 K-1
)
H2
-1.1755
0.0125
0.0014
-1.1620
0.0003
31.135
-1.1768
-0.0563d
-0.4192e
e-(aq)
H+(aq)
-0.0011
-0.0011
a We have used 1 au ) 27.2116 eV and 1 eV ) 23.0605 kcal mol-1 b At 298.15 K (25 °C). c H ) E + kT for gaseous species; H ) E for
.
solvated species. The thermal energy of an electron is (3/2)kT. d The free energy of hydration is -0.0566 au (from ref 38b). e This free energy is
the sum of the free energy of hydration (-0.4181 au) from ref 38a and the Debye-Hu¨ckel energy of the charge-dielectric interaction.
TABLE 2: Molecular Characteristics for the Three Species PC+, PC, and PCH+ a
ONIOM energy (au)b
electrostatic
contributionc
total
ONIOM
thermal
Born
Onsager
entropyc
species
PC+ (Ox)
B3LYP/6-31G
energy (au)d
energy (au)e
energy (au)e
(cal mol-1 K-1
)
-3008.5186
-3008.6688
-3009.0785
2.2539
2.3790
2.2287
-3006.2647
-3006.2899
-3006.8498
0.2862
0.2809
0.2965
-0.7474 ( 0.0053
-0.9459 ( 0.0067
-0.7474 ( 0.0053
-0.0054
-0.0058
-0.0043
158.062
149.013
147.389
PC (red) (pH 7.0)
PCH+ (red) (pH 3.8)
a We have used 1 au ) 27.2116 eV and 1 eV ) 23.0605 kcal mol-1 b The values reported have been computed using an ONIOM2 technique
.
with DFT-B3LYP method being used for the active site, and Amber force field for the remainder protein and solvent water molecules. c Calculated
using the effective dielectric constant ꢀp ) 8.0. d The values reported have been computed using the DFT-B3LYP method for the optimized active
site. e The Born and Onsager energies have been calculated by taking the average radius 42.27 ( 0.30 au for PC and dielectric constant 78.5 for
water at 25 °C.
TABLE 3: Calculation of ∆H° a
pH
reaction
∆EONIOM (au)
∆EThermal (au)
∆EMedium (au)
∆H° (au)a
7.0
3.8
PC+(aq) + e-(aq) f PC(aq)
-0.0252
-0.5851
-0.0067
+0.0089
-0.1990 ( 0.0014b
-0.2308 ( 0.0014
PC+(aq) + H+(aq) + e-(aq) f PCH+(aq)
0.0011b
-0.1460c
a We have used 1 au ) 27.2116 eV and 1 eV ) 23.0605 kcal mol-1
.
b ∆H° ) ∆EONIOM + ∆EThermal + ∆EMedium
.
c The ∆H term is negligibly
small at pH 3.8 and pH 7. d After considering ∆H° for H+(aq), -0.4290 au; ∆H° ) ∆G° + T∆S°, where ∆G° ) -0.4181 au (from ref 38) and
∆S° is -23.01 cal mol-1 K-1 (from ref 46a).
addition of an electron to the Cu atom, accompanied by the
dissociation of the imidazole side chain of the His87 residue
and the protonation of the δ nitrogen on the same residue,
bearing on the determination of reduction potential, as both
electrodes have water as the medium. The only residual medium
interaction term arises form the Debye-Hu¨ckel contribution,
because the ionic strengths of the solution at the two electrodes
are different. This is evaluated to be -0.0011 au. The Debye-
Hu¨ckel contribution toward the free energy of the proton at pH
3.8 is similarly evaluated to be -0.0011 au. The free energy
calculated for process 10 is -0.1695 au (-4.612 eV). See Table
1.
PC+(aq) + H+(aq) + e-(aq) f PCH+(aq)
(14)
For a specific process involving the larger species, ∆G° per
reaction is written as
∆G° ) ∆ESCRF + ∆EBorn + ∆Ethermal + kBT∆n - T∆S (15)
Model Complex. The computed ONIOM energies and other
characteristics of the model complexes are listed in Table 2.
The thermal energies and entropies have been evaluated only
for the DFT optimized active sites as discussed earlier.
The enthalpy changes associated with reactions 13 and 14
are listed in Table 3. As the species remain in the solvent
medium, the PV contribution to enthalpy is exceedingly small
and hence the P∆V term has been neglected. The contribution
of the electronic energy changes to the reduction potential has
been considered to arise only from the changes in the active
site.
where ∆n is again the change in the number of species in gas
phase. The kinetic energy of the aquated electron was taken as
that of the free (gaseous) electron. The standard potential for
the reduction reaction is written as
1
2
E° (in V) ) -∆G° + ∆G° H+(aq)/ H2(g) (in eV) (16)
red
[
(
)
]
When protonation (or deprotonation) is involved, ∆G and the
midpoint potential E will be pH-dependent. Because the
reduction of plastocyanin at pH 3.8 involves the addition of a
proton, the midpoint potential will be given by
Entropy Changes. The change in entropy stems from a
variety of sources such as the change in the protein conformation
and the reorganization of the solvent molecules around the
protein. Hence, the entropy values computed for the isolated
active site can account for only one part of ∆S. It is notoriously
difficult to compute the entropy of a reaction, especially when
one is considering a large molecule. The experimental value
for the entropy change for the net reaction, that is, the reduction
of plastocyanin and the simultaneous oxidation of H2, was
determined by Gray et al.45 The observed value is -103.3 cal
mol-1 K-1. An estimate of the entropy of hydration of a proton
has been obtained from the work of Noyes, who experimentally
Ered (in V) ) E° (in V) - 2.303RTpH (in eV) (17)
red
4. Results and Discussion
Aquated Proton. The free energy of hydration of the proton
and that of the electron have been estimated by Zhan et al. to
be -262.4 and -35.5 kcal/mol, respectively, from a first
principles approach.38a,b The determination of the reduction
potential of the proton has been evaluated using the above
hydration energy. The electronic and thermal energies of gaseous
hydrogen have been evaluated using the 6-31G(d) basis set.
However, the hydration energy of the electron has no direct
determined the value to be -25.36 cal mol-1K-1 46a
the ∆S value for the reduction of plastocyanin was calculated
.
From this,

8014 J. Phys. Chem. B, Vol. 108, No. 23, 2004
Datta et al.
TABLE 4: Calculation of the Reduction Potential of Plastocyanina
midpoint E (mV)
pH
reaction
∆H° (au)
-T∆S (au)
calc ∆G (au)
calc
obs
0.0
7.0
3.8
H+(aq) + e-(aq) f 1/2H2(g)
-0.1129b
PC+(aq) + e-(aq) f PC(aq)
-0.2308 ( 0.0014
-0.1460
+0.0475c
-0.0478f
-0.1267 ( 0.0014d
-0.1289g
376 ( 38
379e
436h
PC+(aq) + H+(aq) + e-(aq) f PCH+(aq)
a We have used 1 au ) 27.2116 eV and 1 eV ) 23.0605 kcal mol-1 b Using free energy of proton hydration and electron hydration from ref
.
43a,b. The Debye-Hu¨ckel correction has been introduced. c Observed value, ref 48. Only +0.0043 au can be accounted for by our calculation on
the active site. Hence, 0.0432 au is due to the reorganization of the protein and surrounding water molecules during the reduction. d Using free
energy of electron hydration -0.0566 au (ref 38b) without the Debye-Hu¨ckel term at pH 7. e Reference 9. f Estimated value that would give rise
to the observed midpoint potential of 436 mV. Out of this value, 0.0051 au can be accounted for from the calculated entropies, -0.0109 au is due
to the change of entropy of the proton on being attached, and -0.0420 au must arise from the reorganization of the protein and the water molecules
surrounding it. g Using the free energy of hydration of electron -0.0566 au. (ref 38b) without the ∆H correction at pH 3.8, and 2.303RT pH )
0.0083 au. h Reference 6c. The value contains -2.303RT pH.
(-100.0 cal mol-1 K-1), because the experimental value of
entropy change for the oxidation of H2 to an aqueous proton is
known (-3.3 cal mol-1 K-1).46b The calculated entropies of
mechanism considered for the reaction above pH 5.5 is
consistent with the reduction potential being independent of pH.
The mechanism considered for the reduction below pH 5.5
agrees with the fact that the reduction potential increases by
about 50-60 mV for each unit drop of pH.
PC+ and PC yields a ∆S value of only -9.05 cal mol-1 K-1
.
Hence the remaining -90.95 cal mol-1 K-1 must account for
the reorganization of water molecules around the active site and
the small change in the protein conformation during the process
of reduction. The entropy decreases as more molecules become
associated on reduction.
The calculated potential (376 ( 38 mV) at pH 7 is in excellent
agreement with the observed one (379 mV). The interaction
with medium and the entropy change are found to play the
crucial role in determining the potential. A similar calculation
leads us to predict the entropy of reduction of PC at pH 3.8 to
be +100.60 cal mol-1 K-1. Thus, the reduction at low pH is
associated with a large increase in entropy. This can be attributed
to the following factors: (i) release of free water molecules on
the attachment of a proton to the imidazole moiety of the
complex; (ii) large change in the dipole moment of the
considered species, which results in a change in water organiza-
tion around the active site; (iii) change in the existing hydrogen
bonding network around the rest of the protein. Our results are
in conformity with the suggestion that the thermodynamic
reasons for the acid transition (protonation in reduced form at
low pH) are entropic for plastocyanins.49
Reduction Potential. The standard free energy changes for
the reduction reaction is written as
∆G° ) H°(PC) - H°(PC+) - T∆Sexp
(18)
red
at pH 7, where ∆Sexp is the experimental value for the reduction
at pH 7. The reduction potential has been calculated for reaction
13. This calculation is shown in Table 4. The calculated potential
(376 ( 38 mV) at pH 7 is in excellent agreement with the
observed one (379 mV).2
Below pH 5.5, the free energy change is given by
∆Gred ) H°(PCH+) - H°(PC+) - G°(H+(aq)) -
T∆S + 2.303RTpH (19)
Despite our apparent success in calculating the redox potential
at pH 7, the present work must be understood as being subject
to several limitations. The latter are as follows:
where the entropy change ∆S remains unknown. Because the
midpoint potential is known (436 mV) at pH 3.8,6 a similar
calculation leads us to predict the entropy of reduction of
plastocyanin at this pH to be +100.60 cal mol-1 K-1. Out of
this, -10.67 cal mol-1 K-1 can be accounted for from the
calculated entropies of the involved species PCH+ and PC
(Table 2), and +23.01 cal mol-1 K-1 is attributable to the
change of entropy of an aquated proton46a on being attached to
the protein moiety. Therefore, 88.26 cal mol-1 K-1 accounts
for the reorganization of water molecules around the active site
and the change in protein conformation during the process of
reduction.
(i) The Debye-Hu¨ckel energy has been calculated for the
aquated proton with unit activity. The linearization of the
Poisson-Boltzmann equation is strictly valid at a much lower
concentration that corresponds to a much shorter reciprocal
length, which makes the use of the Debye-Hu¨ckel expressions
(6) and (7) somewhat dubious. Still, the large dielectric constant
of water gives rise to a Debye-Hu¨ckel contribution that is much
smaller than the other solvation energy contributions. Neverthe-
less, the neglect of the Debye-Hu¨ckel term would cause the
calculated reduction potential to decrease by approximately 60
mV.
(ii) The estimation of the entropy of reduction at pH 3.8 rests
upon the assumption that the structure of the oxidized form of
plastocyanin remains unchanged from the structure found at pH
6.
(iii) The geometry optimizations have been undertaken to an
accuracy of 10-5 au in energy. The subsequent thermal energy
calculation would at most introduce an inaccuracy of only a
few millivolts into the computations.
5. Conclusions
A quantum chemical QM/MM calculation of the reduction
potential of the protein plastocyanin from poplar leaves has been
presented here. An ONIOM methodology has been adopted to
divide the protein into smaller parts, based on their relative
contributions to the reduction potential. Density functional
methodology has been used for the active site, where there are
major changes in electronic structure, whereas the remainder
of the protein has been simulated by UFF and Amber force
fields during the optimization of geometry and the single point
energy evaluation, respectively. The fundamental structural
features of the copper-containing protein have thus been retained
during the course of geometry optimization of the species. The
Additionally, reduction potential calculations are found to be
sensitive to the estimated radius of the protein. Thus, a change
of (0.1 au in the protein radius affects the change in Born
energy by about (0.0005 au. This causes the computed
reduction potential to shift by (13 mV. The calculation of the
entropy of reduction of the protein at pH 3.8 is free from this
error, because the net contribution of the Born term is zero in

Standard Reduction Potential of Plastocyanin
J. Phys. Chem. B, Vol. 108, No. 23, 2004 8015
form of the protein. Whereas the precise determination of the
dielectric constant of a protein is still not possible, this method
allows for an estimation of the effective dielectric constant of
the protein in solution. Simonson et al. have discussed that there
is an error of about 3 in the estimation of dielectric constant by
a molecular dynamics technique.
We undertook a 100 ps simulation on a section of the protein
to estimate the approximate degree of penetration of water into
the protein. Altogether 640 water molecules were taken into
consideration. The MM+ force field was chosen for the
molecular dynamics simulation, as it is an all atom force field,
which enables placing of quantum mechanically computed
charges on the atoms. This force field includes cutoffs for
nonbonded interactions, which makes it a very general force
field for solvation and molecular dynamics simulations.47 The
protein was first optimized using an MM+ force field and the
DFT generated charges on the solvent (water) molecules and
the protein. The MM+ force field was then used to simulate
the dynamics of a section of the protein in water at 300 K, in
steps of 0.002 ps. Hyperchem Professional Release 7 for
Windows47 was used in these calculations. The results show
that the water molecules effectively penetrate into a layer of
thickness in range of 4-10 au into the protein. This is shown
by the shaded region in Figure 6.
Figure 6. Schematic diagram for the protein in water, following ref
43a,b.
this case. Because of these reasons, the consideration of the
Born energy as a first order correction suffices.
A more important factor is the dielectric constant in the
interior of the solvated protein. The value of 8 adopted in this
work has led to a more or less correct reduction potential at pH
7.0. This is rather fortuitous. Nevertheless, if ꢀp changes by one
unit from the value 8, the calculated potential changes by about
400 mV. Therefore, a very high value of ꢀp in the bulk of
plastocyanin is entirely unrealistic.
Simonson et al. have discussed that the outer layer generally
has an effective dielectric constant of 20. The dielectric constant
of dry cytochrome c powder is 3.6.44b Hence the average
dielectric constant can be estimated as
In short, the calculated electronic and solvent characteristics
are in good agreement with the electronic structure of the Cu-
containing region in the real specimen, plastocyanin. Determi-
nation of the rate of electron transfers from cytochromef to
plastocyanin and from plastocyanin to P700+ is our ultimate
objective and this work prepares the basis for our future work
in this direction.
3
3
r1 ꢀ1 + (rg3 - r1 )ꢀ2
ꢀp )
(A.1)
3
rg
Acknowledgment. We thank the reviewers for making a
large number of critical and valuable suggestions.
As discussed in the text, the effective radius of the protein (rg)
is 42.27 au. Therefore, the average dielectric constant in the
bulk is found in the region of 7.8-12.7 for rg - r1 ) 4-10 au.
The value of 12.7 is certainly an overestimation, as only a few
water molecules penetrate up to a thickness of 10 au. The lower
average value 8 has been used in the text, as most of the solvent
penetration extends up to approximately 4 au.
Appendix I: Evaluation of the Dielectric Constant of the
Interior of a Protein
The explicit derivation of an effective dielectric constant (ꢀ1)
for the inner side of a globular protein is described here.
Molecular dynamics has provided several insights into the effect
of charge screening in proteins.43a,b Results from molecular
dynamics indicate that the dielectric constant of a protein can
be hardly understood from the dielectric constant values obtained
from the physical measurements undertaken on dry powders of
proteins. Two important factors that affect the dielectric constant
of a protein in solution are (i) the ionization of residues
containing carboxylate and amino groups and (ii) the interaction
of solvent with the outer layers of the protein. Following the
molecular dynamical investigations in ref 43a,b, a globular
protein can be modeled as a sphere with two concentric shells
having different dielectric constants. For a water soluble protein
such as plastocyanin, the outer layer is highly polar, with
dielectric constant ꢀ2. We take ꢀp as the net dielectric constant
of the protein, which is to be determined. Besides, rg is the
effective radius of the whole entity, and r1 is the effective radius
of the inner shell, as shown in Figure 6.
As shown in Figure 2, the active site is located in the outer
layer. The charged residues are also located in the outer
periphery of the protein. Consequently, the effective dielectric
constant of the protein must be determined and considered for
the evaluation of the energy of interaction of the charge and
the dipoles in the bulk with the active site. The effect and the
degree of solvation of the outer peripheral layer of the protein
were estimated from a molecular dynamics study of the oxidized
References and Notes
(1) Redinbo, M. R.; Cascio, D.; Chonkair, M. K.; Rice, D.; Merchant,
S.; Yeates, T. O. Biochemistry 1993, 32, 10560.
(2) Katoh, S.; Shiratori, I.; Tahamiya, A. J. Biochem. (Tokyo) 1962,
51, 32.
(3) Fee, A. J. Structure Bonding; Springer: Berlin, 1975; Vol. 23, p
60. Colman, P. M.; Freeman, H. C.; Guss, J. M.; Murata, M.; Norris, V.
A.; Ramshaw, J. A. M.; Venkatappa, M. P. Nature 1978, 272, 319.
(4) Sykes, A. G. AdV. Inorg. Chem. 1991, 36, 377. Machonkin, T. E.;
Zhang, H. H.; Hedman, B.; Hodgson, K. O.; Solomon, E. I. Biochemitsry
1998, 37, 9570.
(5) Guss, J. M.; Harrowell, P. R.; Murata, M.; Norris, V. A.; Freeman,
H. C. J. Mol. Biol. 1986, 192, 361.
(6) (a) Beinert, H. Coord. Chem. ReV. 1980, 33, 55. (b) Guss, J. M.;
Freeman, H. C. J. Mol. Biol. 1983, 169, 521. (c) Buchi, F. N.; Bond, A.
M.; Codd, R.; Huq, L. N.; Freeman, H. C. Inorg. Chem. 1992, 31, 5007.
(7) Segal M. G.; Sykes, A. G. J. Am. Chem. Soc. 1978, 100, 4584.
(8) (a) Markley, J. L.; Ulrich, E. L.; Krogmann, D. W. Biochem.
Biophys. Res. Commun. 1977, 78, 106. (b) Freeman, H. C.; Norris, V. A.;
Ramshaw, J. A. M.; Wright, P. E. FEBS Lett. 1978, 86, 131.
(9) Datta, S. N.; Priyadarshy, S. Chem. Phys. Lett. 1990, 173, 360.
(10) Mallik, B.; Datta, S. N. Int. J. Quantum Chem. 1994, 52, 629.
(11) Datta, S. N.; Mallik, B. Int. J. Quantum Chem. 1995, 53, 37; 1997,
61, 865.
(12) Datta, S. N.; Mallik, B. J. Phys. Chem. B 1997, 101, 4171, 5191.
Datta, S. N.; Prabhakar, B. G. S.; Nehra, V. J. Phys. Chem. B 1999, 103,
2291. Datta, S. N.; Nehra, V.; Jha, A. J. Phys. Chem. B 1999, 103, 8768.
Datta, S. N.; Parandekar, P. V.; Lochan, R. C. J. Phys. Chem. B 2001, 105,
1442.

8016 J. Phys. Chem. B, Vol. 108, No. 23, 2004
Datta et al.
(13) Datta, S. N.; Deshpande, R. Ind. J. Pure Appl. Phys. 1997, 35, 483.
(14) (a) Garrett, T. P. J.; Clingeleffer, D. J.; Guss, J. M.; Rogers, S. J.;
Freeman, H. C. J. Biol. Chem. 1984, 259, 2822. (b) Guckert, J. A.; Lowery,
M. D.; Solomon, E. I. J. Am. Chem. Soc. 1995, 117, 2817. (c) De Kerpel,
Jan, O. A.; Ryde, U. Proteins: Struct., Funct., Genet. 1999, 36 (2), 157.
(15) Kurnikov, I. V. Beratan, D. N. J. Chem. Phys. 1996, 105, 9561.
(16) Newton, M. D. Int. J. Quantum. Chem. Symp. 1980, 14, 363.
(17) Logan, J.; Newton, M. D. J. Chem. Phys. 1983, 78, 4086.
(18) Balaji, V.; Jordan, K. D.; Burrow, P. D.; Paddon-Row, M. N.;
Patney, H. K. J. Am. Chem. Soc. 1982, 104, 6849.
(19) Jordan, K. D.; Paddon-Row, M. N. Chem. ReV. 1992, 92, 395.
(20) Mikkelsen, K. V.; Ratner, M. A. J. Phys. Chem. 1989, 93, 1759.
(21) Gruschus, J. M.; Kuki, A. J. Phys. Chem. 1993, 97, 5581.
(22) Larsson, S.; Volosov, A. J. Chem. Phys. 1987, 87, 6623.
(23) Beratan, D. N.; Hopfield, J. J. J. Am. Chem. Soc. 1984, 106, 1584.
(24) Beratan, D. N.; Onuchic, J. N.; Hopfield, J. J. J. Chem. Phys. 1987,
86, 4488.
Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A.
D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi,
M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.;
Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick,
D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.;
Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi,
I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.;
Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M.
W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L.; Head-Gordon,
M.; Replogle, E. S.; Pople, J. A. Gaussian 98, revision A.11.4; Gaussian,
Inc.: Pittsburgh, PA, 1998.
(37) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K.
M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.;
Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 8897.
(38) (a) Zhan, C. G.; Dixon, D. A. J. Phys. Chem. A 2001, 105, 11534.
(b) Zhan, C. G.; Dixon, D. A. J. Phys. Chem. B 2003, 107, 4403.
(39) Newton, M. D. J. Phys. Chem. 1975, 79, 2795.
(40) Olsson, M. H. M.; Ryde, U. J. Biol. Inorg. Chem. 1999, 4, 654.
(41) Tomasi, J.; Persico, M. Chem. ReV. 1994, 94, 2027.
(42) (a) Torres, R. A.; Lovell, T.; Noodleman, L.; Case, D. A. J. Am.
Chem. Soc. 2003, 125 (7), 1923. (b) Olsson, M. H. M.; Hong, G.; Warshel,
A. J. Am. Chem. Soc. 2003, 125, 5025. (c) Johnson, E. T.; Parson, W. W.
Biochemistry 2002, 41 (20), 6483. (d) Claudia, S. N.; Warshel, A.
Proteins: Struct., Funct., Genet. 2001, 44 (4), 400.
(43) (a) Simonson, T.; Brooks, C. L., III. J. Am. Chem. Soc. 1996, 118,
8452. (b) Simonson, T.; Perahia, D. Proc. Natl. Acad. Sci. U.S.A. 1995,
92, 1082.
(44) Bockris, J. O’M; Reddy, A. K. N. In Modern Electrochemisrty;
Plenum Press: New York and London, 1998; p 565.
(45) Napanon, S.; Anson, F. C.; Gray, H. B. J. Am. Chem. Soc. 1979,
101, 455.
(25) Priyadarshy, S.; Skourtis, S. S.; Risser, S. M.; Beratan, D. N. J.
Chem. Phys. 1996, 104, 9473.
(26) Larsson, S. J. Am. Chem. Soc. 1981, 103, 4034; J. Chem. Soc.,
Faraday Trans. 2 1983, 79, 1375; Chem. Phys. Lett 1982, 90, 136; J. Phys.
Chem. 1984, 88, 1321.
(27) Ratner, M. A. J. Phys. Chem. 1990, 94, 4877.
(28) Broo, A.; Larsson, S. J. Phys. Chem. 1991, 95, 4925.
(29) Stuchenbrukhov, A. A.; Marcus, R. A. J. Phys. Chem. 1995, 99, 7581.
(30) Evenson, J. W.; Karplus, M. J. Chem. Phys. 1992, 96, 5272.
(31) Ryan, J. J.; Risser, S. M.; Beratan, D. N.; Onuchic, J. N. J. Phys.
Chem. 1993, 97, 13083.
(32) Siddarth, P.; Marcus, R. A. J. Phys. Chem. 1992, 96, 3213; J. Phys.
Chem. 1993, 97, 2400; J. Phys. Chem. 1993, 97, 6111.
(33) Svensson, M.; Humbel, S.; Froese, R. D. J.; Matsubara, T. J. Phys.
Chem. 1996, 100, 19357.
(46) (a) Noyes, R. M. J. Am. Chem. Soc. 1962, 84 (4), 513. (b)
Battistuzzi, G.; Borsari, M.; Canters, G. W.; Waal, E. de; Leonardi, A.;
Ranieri, A.; Sola, M. Biochemistry 2002, 41, 14293.
(47) HyperChem Professional Release 7.01 for Windows; Hypercube
Inc., Gainesville, FL, 2002.
(34) Rappe, A. K.; Casewit, C. J.; Colwell K. S.; Goddard, W. A., III;
Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024.
(35) Rappe, A. K.; Goddard, W. A., III. J. Phys. Chem. 1991, 95, 3358.
(36) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb,
M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A., Jr.;