Acidity of a Cu-Bound Histidine in the Binuclear Center of Cytochrome c Oxidase
Elisa Fadda, Nilmadhab Chakrabarti, and Re´gis Pome`s*
Structural Biology and Biochemistry, The Hospital for Sick Children, and Department of Biochemistry,
UniVersity of Toronto, Toronto, Ontario, Canada
ReceiVed: May 24, 2005; In Final Form: August 16, 2005
Cytochrome c oxidase (CcO) is a crucial enzyme in the respiratory chain. Its function is to couple the reduction
of molecular oxygen, which takes place in the Fea3-CuB binuclear center, to proton translocation across the
mitochondrial membrane. Although several high-resolution structures of the enzyme are known, the molecular
basis of proton pumping activation and its mechanism remain to be elucidated. We examine a recently proposed
scheme (J. Am. Chem. Soc. 2004, 126, 1858; FEBS Lett. 2004, 566, 126) that involves the deprotonation of
the CuB-bound imidazole ring of a histidine (H291 in mammalian CcO) as a key element in the proton pumping
mechanism. The central feature of that proposed mechanism is that the pKa values of the imidazole vary
significantly depending on the redox state of the metals in the binuclear center. We use density functional
theory in combination with continuum electrostatics to calculate the pKa values, successively in bulk water
and within the protein, of the Cu-bound imidazole in various Cu- and Cu-Fe complexes. From pKas in bulk
water, we derived a value of -266.34 kcalâmol-1 for the proton solvation free energy (¢GsolH
+). This estimate
is in close agreement with the experimental value of -264.61 kcalâmol-1 (J. Am. Chem. Soc. 2001, 123,
7314), which reinforces the conclusion that ¢GsolH
+ is more negative than previous values used for pKa
calculations. Our approach, on the basis of the study of increasingly more detailed models of the CcO binuclear
center at different stages of the catalysis, allows us to examine successively the effect of each of the two
metals’ redox states and of solvation on the acidity of imidazole, whose pKa is approximately 14 in bulk
water. This analysis leads to the following conclusions: first, the effect of Cu ligation on the imidazole
acidity is negligible regardless of the redox state of the metal. Second, results obtained for Cu-Fe complexes
in bulk water indicate that Cu-bound imidazole pKa values lie within the range of 14.8-16.6 throughout
binuclear redox states corresponding to the catalytic cycle, demonstrating that the effect of the Fe oxidation
states is also negligible. Finally, the low-dielectric CcO proteic environment shifts the acid-base equilibrium
toward a neutral imidazole, further increasing the corresponding pKa values. These results are inconsistent
with the proposed role of the Cu-bound histidine as a key element in the pumping mechanism. Limitations
of continuum solvation models in pKa calculations are discussed.
1. Introduction
Cytochrome c oxidase (CcO) is the terminal enzymatic
complex in the mitochondrial respiratory chain. Its catalytic
function is to reduce O2 to H2O and to harness the free energy
released by that reaction to pump protons across the mitochon-
drial inner membrane. The electrochemical gradient generated
across the membrane by the excess protons is ultimately used
to drive ATP synthesis.1-5 The high efficiency of CcO rests
precisely on the balanced coupling between proton pumping
and the redox reaction:
The oxygen reduction takes place in the binuclear center, which
consists of a high-spin Fe hemea3 and a CuB, in the enzyme
active site (see Figure 1). Two distinct proton uptake pathways
connect the binuclear center to the matrix side of the mitochon-
drial membrane,5-8 namely the D-channel and the K-channel,
which owe their names to conserved aspartate and lysine
residues, respectively. It is generally accepted that seven of the
eight protons involved in the catalytic cycle follow the D-
channel from the matrix to the active site, while only one proton
adopts the K-channel.9,10 Four of the eight protons are consumed
in the redox reaction to produce two water molecules (chemical
protons), while the other four are translocated through the
enzyme (pumped protons). All protons are transported to the
CcO active site, from which they are directed either toward the
binuclear center or toward the exit pathway. A detailed picture
of the mechanism underlying this selective proton distribution
is a matter of debate.11-15 The analysis of the many available
crystallographic structures of CcO does not show any specific
pathway connecting the mouth of the D-channel to the binuclear
center or leading protons to the exterior side of the membrane.
It is thought that both chemical and pumped protons are relayed
within the active site cavity through dynamic networks of
hydrogen-bonded water molecules.7,13,16-18 Furthermore, the
direction of proton transfer through hydrogen-bond networks
has been suggested to be coupled to the electron transfer to the
hemea3-CuB binuclear center.13 Recent studies have proposed
the direct involvement of one of the CuB histidine ligands (i.e.,
H291 in bovine heart CcO, or H334 in Rhodobacter sphaeroides
CcO) as a proton loading site in the pumping mecha-
nism.7,11,14,15,19 In these mechanisms, the protonation state of
Nä1 in H291 is assumed to depend on the redox state of the* Corresponding author. E-mail: pomes@sickkids.ca.
O2 + 4e
- + 8H+ f 2H2O + 4H
+ (1)
22629J. Phys. Chem. B 2005, 109, 22629-22640
10.1021/jp052734+ CCC: $30.25 © 2005 American Chemical Society
Published on Web 11/08/2005

metals in the binuclear center. Wikstro¨m proposed that the
protonation and deprotonation of the histidine side chain is
coupled to the dissociation of its imidazole ring from CuB.11
An alternative scheme7,14,15,19 entails that H291 remains ligated
to CuB and that deprotonation at the Nä1 position causes the
formation of an imidazolate anion (i.e., CuB-ImH f CuB-
Im-). This hypothesis was proposed on the basis of continuum
electrostatic calculations of CcO.14,15 A key feature of that
proposal is that the entry of a chemical proton in the active
site, which converts an OH- ligand of CuB to H2O, causes the
pKa of H291 to drop, leading to the deprotonation of Nä1.
Although the deprotonation of imidazole, which produces an
imidazolate anion, is not possible in bulk water (pKa  14),20-22
the presence of a metal-coordinated imidazolate anion has
recently been determined experimentally in the iron-sulfur
protein of the cytochrome bc1 complex.23 In this system, the
imidazolate ring of one of the two Fe-coordinated histidines is
hydrogen bonded to the ubiquinol substrate and presumably
represents the redox-linked protonation site. A transient endiol-
imidazolate complex was also identified as a probable inter-
mediate in the catalytic mechanism of triosephosphate isomerase.24
The likelihood of the formation of a Cu-bound imidazolate in
the CcO active site can be addressed computationally by using
a suitable quantum mechanical (QM) approach on an appropriate
model system. In this study, we use density functional theory
(DFT) in combination with continuum electrostatics to examine
the dependence of the Cu-bound imidazole pKa on the redox
state of different Cu and Cu-Fe complexes. Our approach is
based on the analysis of the Cu-bound imidazole acidity in
increasingly larger systems modeled after the CcO binuclear
center at various stages of the catalytic cycle. A schematic
depiction of the main steps of the reduction of O2 in the
binuclear center, and the proposed connection to the H291
deprotonations,14 are shown in Figure 2. The fully reduced state
of the binuclear center, defined as the “R” state, will be
considered here as the starting point of the catalytic cycle. An
O2 enters the active site and binds to the metals, forming a
“peroxy” (P) state (not shown). The tyrosine (Y244) cross-linked
to His240 forms a radical by releasing a hydrogen atom, which
helps cleaving the peroxy bond.25-27 The first electron transfer
from the hemea allows the formation of a negatively charged
tyrosine (indicated in Figure 2 as YO- in the Pr state). In the
“ferryl” state (F) the first chemical proton (HIn+) accesses the
binuclear center and saturates Y244. The latter remains neutral
throughout the following steps of the catalytic cycle.28 In the
transition between the F and Fp states, a chemical proton enters
the binuclear center and converts the hydroxyl CuB ligand into
water, increasing the total charge of the site from +1 to +2.
That stage of the catalytic cycle is supposedly connected to the
first deprotonation of the His291 imidazole ring. The next proton
entering the active site reloads the His291 Nä1 position. A
dismutation leads to the formation of the H′ state. The transition
between the H′ and O states is connected to a second
deprotonation of H291. In the E state, the first water molecule
produced in the redox reaction leaves the binuclear center and
the redox state of Fea3 changes from III to II. Here, the fourth
and last chemical proton enters the active site and saturates the
hydroxyl (Ep state), triggering the third His291 deprotonation.
The next incoming proton reloads the Nä1 position, closing the
catalytic cycle. In the original article proposing the deprotonation
of His291, the catalytic step connected to the fourth pumping
event was not clearly identified.14
The complexes studied, shown schematically in Figure 3,
were designed to represent stages of the CcO catalytic cycle
presumably connected to the Cu-bound imidazole deprotona-
tions. First, we determined the pKa values in bulk water for
simple [ImH]3CuM-L models (where L indicates the fourth
ligand, i.e H2O or OH-, and M indicates the oxidation state of
Cu, i.e., I or II). The results on these complexes show the effect
of the redox states of the central Cu on the acidity of an isolated
imidazole. The effect of the Fe porphyrin was added succes-
sively by studying more complex systems where various Fe
oxidation states (i.e. II, III, IV) and various intermetal ligands
were considered. Last, we calculated the pKa of the imidazole
ring as a CuB ligand in the binuclear center of CcO to assess
the effect of the low dielectric proteic environment.
The calculation of absolute pKa values of species in solution
is a delicate task, which necessitates the use of high-level QM
methods capable of attaining chemical accuracy (average
absolute error of about 0.1 eV (or 2 kcalâmol-1) against the G2
thermochemical data set.29). The treatment of solvation repre-
sents the least reliable aspect of pKa calculations. The expected
limits of accuracy range within 2-5 kcalâmol-1.30 To gauge
Figure 1. Subunit A of cytochrome c oxidase (CcO) from Rhodobacter sphaeroides (PDB structure 1M56, ref 6). The residues in the binuclear
center are shown in detail at the top-right corner. The numbering of all residues correspond to mammalian CcO; Cu is shown in teal and Fe in
orange. The location of the active site within CcO is indicated by an arrow.
22630 J. Phys. Chem. B, Vol. 109, No. 47, 2005 Fadda et al.

the dependence of pKa values on the level of theory used to
describe solvation, an SCF generalized Born (GB) model and
a post-SCF Poisson-Boltzmann (PB) approach are used to
determine the solvation free energy in bulk water. Furthermore,
tests on acid-base couples for which the solvation free energy
is known experimentally show the limitations of GB and PB
methods in the treatment of ionic species.
The Cu-bound histidine pKa values are calculated within CcO
by using PB continuum electrostatic calculations. In this highly
inhomogeneous environment, the pKas must be determined
through a thermodynamic scheme that includes explicitly the
contribution of proton solvation free energy (¢GsolH
+). Many
experimental and computational studies31-35 have focused on
the determination of ¢Gsol
H+
, and the values reported therein
cover a range of 12.0 kcalâmol-1 (from -252.6 to -264.6
kcalâmol-1). Because of this large range of uncertainty, and for
consistency in our results, we estimated ¢Gsol
H+ based on the
pKa values that we calculated in bulk water (see Section 2 for
details). We obtained ¢GsolH
+
) -266.340 kcalâmol-1, close to
the value of -264.61 kcalâmol-1 published by Liptak and
Shields,31 which corresponds to the lowest value in the range
of the proton solvation free energies reported in the literature.
Our estimate reinforces the conclusion that ¢Gsol
H+ is more
negative than previous values used for pKa calculations.
Our calculations show that the acidity of the CuB-bound
imidazole ring in the complexes studied does not depend
significantly on the redox state of the metals. The pKas of the
Cu-bound imidazole ring in bulk water vary between 14.8 and
16.6 for the large Cu-Fe complexes around the catalytic cycle.
The low dielectric environment of the protein shifts the
equilibrium toward the protonated (neutral) form of the Cu-
bound histidine, increasing the pKa values. The implications of
these results on the proposed CcO pumping mechanism are
analyzed. In two recent publications,36,37 the pKa values of
imidazole ligands in Cu complexes and in CuB-hemea3 models
were evaluated by using an approach similar to ours. However,
fundamental differences in the pKa calculation methodology,
which will be discussed in detail in the following sections, lead
us to radically different conclusions.
Figure 2. Proposed main steps of the reduction of O2 to H2O in the binuclear center of CcO. The hemea3 is represented by the central Fe atom and
by a line defining the porphyrin ring. Only the CuB histidine ligand side chain undergoing deprotonation is shown. The protonated and deprotonated
forms of imidazole are indicated as ImH and Im-, respectively. Catalytic states are defined by a letter shown on top of each frame, according to
the commonly used nomenclature.11,15 The protons entering and exiting the binuclear center are indicated as HIn
+
and HOut
+
, respectively.
Figure 3. Left Cu- and right Cu-Fe complexes studied. The catalytic
states of CcO corresponding to the Cu-Fe models are indicated at the
top-left corner.
Acidity of a Cu-Bound Histidine in Cytochrome c Oxidase J. Phys. Chem. B, Vol. 109, No. 47, 2005 22631

2. Computational Method
2.1. pKa Calculations from First Principles. The thermo-
dynamic cycle commonly used to calculate the acidity of a
chemical group AH is shown in Figure 4. On the basis of that
cycle, the pKa can be expressed through the equation:
where T ) 298 K and R is the molar gas constant. The gas-
phase deprotonation free energy (¢GgasAH) is given by
It comprises the gas-phase deprotonation energy (¢HgasAH), a
vibrational term (¢GThA H), the proton translational energy
(HtransH
+
) 3/2kBT), a term accounting for the change of volume
during the reaction (¢(pV) ) kBT), and the entropic contribution
(SH+) to the proton gas-phase energy (7.8 kcalâmol-1 according
to the Sackur-Tatrode equation38). ¢HgasAH and ¢GThAH can be
evaluated to high accuracy by first principle calculations. The
solvation free energy (¢¢GsolAH in eq 2) is given by
Both ¢Gsol
A-
and ¢Gsol
AH
are usually evaluated through self-
consistent continuum solvation methods.31,39 The proton sol-
vation free energy (¢GsolH
+) represents an elusive quantity,
which is very difficult to determine experimentally and to
calculate reliably.40,41 Indeed, the available experimental31-34
and calculated35 values reported for ¢Gsol
H+
vary between
-252.6 and -264.6 kcalâmol-1. A deviation of only 1.4
kcalâmol-1 in any of the free energy contributions appearing in
the numerator of eq 2 causes the pKa to shift by 1 unit. Hence,
by taking into account ¢Gsol
H+ in the thermodynamic cycle, one
can only determine acidities with an uncertainty of ap-
proximately 9 pKa units. Ultimately, the pKas determined
through this method are linearly dependent on the choice of
the “best” value for ¢Gsol
H+
.
The use of an alternative scheme that circumvents the problem
of accounting for proton solvation has proven to be quite
successful, leading to an accuracy within half of a pKa unit for
carboxylic acids.42 Within this framework, shown in Figure 5,
the unknown pKa of an acid AH is determined relative to that
of the known pKa of another species, BH
where
and
The proton solvation free energy is implicitly taken into account
by the introduction of pKa
BH
as a known constant in eq 5. The
pKa values calculated according to this scheme are independent
of the chosen reference (BH). Therefore, in principle, any type
of acid BH can be selected. However, limitations in the accuracy
of continuum solvation models in describing ionic species
accurately can be the cause of additional errors.
2.2. Model Systems. The pKa values of an imidazole ligated
to Cu were determined in several complexes modeled after the
binuclear center of CcO in catalytic stages proposed14 to be
connected to H291 deprotonation (see Figure 3). We first
calculated the pKa of imidazole (ImH) as a ligand in [ImH]3CuML-
type complexes, where the superscript M indicates the oxidation
state of the metal, while L indicates the fourth ligand (H2O or
OH-). From now on, we will refer to this type of complex as
the “Cu model”. From the CcO catalytic cycle shown in Figure
2, we selected three Cu models for the pKa calculation,
namely: [ImH]3Cu IH2O, [ImH]3Cu IIH2O (shown as structure
a in Figure 6), and [ImH]3Cu IIOH.
The effect of the Fe porphyrin on the pKa at a close distance
to the Cu center (4.82 Å in Rhodobacter sphaeroides CcO, PDB
code 1M56, ref 6) is then studied in larger model complexes,
which we will refer to as the “Cu-Fe models”. [ImH]3 Cu II-
H2OâââOdFe IV[(porph)NH3] is modeled after the F-state of the
binuclear center. Two systems are modeled on the O-state,
[ImH]3 Cu II-OHâââH2O-Fe III[(porph)NH3] (shown as struc-
ture b in Figure 6), and [ImH]3 Cu II-OHâââFe III[(porph)NH3],
respectively, with and without the water molecule ligated to
Fe. For the latter model, both antiferromagnetic (S ) 2), and
ferromagnetic (S ) 3) coupling between the two metals are
considered to assess the influence of spin coupling on ImH
acidity. Last, [ImH]3 Cu II-OHâââFe II[(porph)NH3] represents
the E state. In all complexes, the apical ligand NH3 of the heme
group replaces a His side chain. For both Cu-Fe and Cu models,
the acidic proton belongs to Nä1 (indicated by an asterisk in
Figure 6). The imidazole that is deprotonated corresponds to
H291 in the CcO binuclear center.
Calculations in bulk water were carried out by using the
relative scheme depicted in Figure 5. The evaluation of the
solvation free energy (¢¢Gsolref in eq 5) represents the least
reliable step in this calculation, especially for the characterization
of ionic species. To assess the dependence of the pKa values
on the solvation model used, we determined ¢¢Gsol
ref for acid-
base couples for which ¢¢Gsol
ref is known experimentally. The
results obtained with two continuum models, a self-consistent
Figure 4. Thermodynamic cycle used to determine the pKa of the
species AH. ¢Gsol represents the free energy of solvation.
Figure 5. Relative scheme for the calculation of the pKa of species
AH. pKa1 is the unknown value to determine, pKa2 is the known
experimental value of the species BH.
¢¢Gsol
ref ) (¢GsolA
-
- ¢Gsol
AH) - (¢GsolB
-
- ¢Gsol
BH) (7)
pKa
AH )
(¢GgasAH + ¢¢GsolAH)
2.303RT (2)
¢Ggas
AH ) ¢Hgas
AH + ¢GTh
AH + Htrans
H+ + ¢(pV) - T[SH+] (3)
¢¢Gsol
AH ) ¢Gsol
A- + ¢Gsol
H+ - ¢Gsol
AH (4)
pKa
AH ) pKa
BH +
(¢¢Ggasref + ¢¢Gsolref)
2.303RT (5)
¢¢Ggas
ref ) (¢HgasAH + ¢GThAH) - (¢HgasBH + ¢GThBH) (6)
22632 J. Phys. Chem. B, Vol. 109, No. 47, 2005 Fadda et al.

generalized Born (GB) model and a Poisson-Boltzmann (PB)
post-SCF approach, were compared to the experimental data.
This analysis allows us to gauge the level of accuracy attainable
for the pKa value in neutral and ionic species. Furthermore, it
allows us to choose as reference equilibrium for the relative
scheme (i.e., BH f B- + H+, see Figure 5) the acid-base
couple whose solvation energy is described most accurately by
both solvation models.
For the calculations in CcO, one cannot use the relative
scheme shown in Figure 5, simply because the pKa of any acid
that could serve as reference BH is not known in the CcO proteic
environment. Therefore, we have to rely instead on the ther-
modymamic cycle depicted in Figure 4. Accordingly, the pKas
are determined through eq 2 and 3. Although the calculations
of the imidazole acidity in the protein interior are designed to
reproduce most closely the effect of such a complex environ-
ment, it is necessary to keep in mind that these pKa values are
linearly dependent on the choice of the “best” value for the
proton solvation free energy. In view of this problem and for
consistency of the results, we chose to determine ¢Gsol
H+
via eq
4 based on the pKas obtained for Cu and Cu-Fe models in
bulk water. More specifically, ¢Gsol
H+
was obtained by plugging
into eq 4 calculated values for ¢Gsol
A-
and ¢Gsol
AH
, together with
the value of ¢¢Gsol
AH
obtained from the thermodynamic cycle
depicted in Figure 5. The value we obtained (-266.340
kcalâmol-1) was then used to evaluate the Cu-bound imidazole
ring acidity in the protein.
2.3. Calculation Details. The initial coordinates for all the
models studied were taken from the X-ray structure of CcO
from Rhodobacter sphaeroides, PDB code 1M56.6 The different
ligands (i.e., H2O, OH-, and O) present at different stages of
the catalytic cycle were added manually to the initial structures.
The coordinates of all Cu-Fe and Cu models are given as
Supporting Information. All open-shell calculations were per-
formed with unrestricted hybrid DFT, namely with the B3LYP
exchange-correlation (xc) functional (U-B3LYP).44 Because of
the small influence of spin contamination that may be introduced
by the exact HF exchange in B3LYP, it is important to underline
that, for all SCF-converged calculations, the expectation value
of the total spin 〈S2〉 differs from S(S + 1) by less than 10%.
All protonated and deprotonated Cu models were fully
optimized at the 6-31G(1p,1d) level. The electrostatic gas-phase
contributions and the vibrational analysis were performed with
an all-electron TZV(1p,1d) basis set. For the Cu-Fe complexes,
all degrees of freedom of the protonated species were relaxed
to a gradient convergence threshold of 0.01. This initial step
allowed the internal structure to adapt to the addition of the
ligands and to the redox state of the metals. This procedure did
not lead to significant deviations from the original X-ray
structure of the binuclear center. For these large complexes, the
deprotonation is assumed to be instantaneous; therefore, the
structure of the deprotonated form is assumed to be identical
to that of the protonated form. This assumption is supported by
the good agreement observed between the pKa values obtained
for the large Cu-Fe complexes and the Cu models, which were
optimized in both protonated and deprotonated forms (see
Section 3).
Geometry optimizations on the Cu-Fe models were per-
formed with SBKJC effective core potentials (ECP) for all heavy
atoms,45,46 while valence electrons were described at the 4-31G
level. The electrostatic gas-phase terms were determined with
a LACV3P(1p,1d) basis set,47 which entails ECPs for Cu and
Fe, and a TZV(1p,1d) quality basis set for all electrons in H,
C, N, and O atoms, and valence electrons in the metals. The
¢GTh terms were assumed to be equal to the ones calculated
for the Cu models. Because deprotonation can be considered
as a relatively localized phenomenon, and because the Cu model
represents a wide enough portion of the Cu-Fe model, this
approximation is reasonable. Solvation free energies for Cu
models were determined by using two different levels of
theory: a self-consistent generalized Born (GB) approximation
based on CM3 model charges,48,49 calculated at the B3LYP//
6-31G(1d) level, and a continuum electrostatic method based
on the Poisson-Boltzmann (PB) equation.50 The latter requires
the calculation of the point charges fitted to the electrostatic
potential (ESP). Here, these point charges are determined with
the Singh-Kollman method,51 also at the B3LYP//6-31G(1d)
level. This protocol was chosen for the sake of consistency with
the level of theory used for the GB-CM3 calculations, to which
the PB results are compared. Extensive testing52 demonstrated
that it yields accurate results compared to that of higher
correlated methods. All ESP charges in the QM region were
calculated in the gas phase. This approach is consistent with
the evaluation of electrostatic parameters for empirical force
fields, most of which are based on ab initio calculations in the
gas phase.53-57 The charge distribution is placed in a cavity
with QM ) 1, which is immersed in a continuum dielectric of
 ) 78.4 representing bulk water. The solvation free energy is
then determined by solving the PB equation by using a finite
difference method. The atomic radii used to construct the solvent
accessible surface (SAS) for H, C, N, and O atoms correspond
to the van der Waals (vdW) radii from Bondi.58 The vdW radius
Figure 6. (a) [ImH]3CuII-H2O complex, and (b) [ImH]3CuII-OHâââH2O-FeIII[(porph)NH3] complex. Cu and Fe atoms are shown in teal and in
orange, respectively. The imidazole which undergoes deprotonation is indicated by an asterisk.
Acidity of a Cu-Bound Histidine in Cytochrome c Oxidase J. Phys. Chem. B, Vol. 109, No. 47, 2005 22633

for Fe was chosen as 1.456 Å.59 Because both metals are buried
inside the complexes, their vdW spheres do not contribute to
the formation of the SAS, so that the precise choice of their
vdW radii is not crucial to the accuracy of the solvation free
energy. The solvation free energy of the Cu-Fe models in bulk
water was determined only with the PB method due to SCF
convergence difficulties in the GB-CM3 calculations.
The deprotonation of the CuB-bound histidine (at pH ) 7)
within CcO was studied by using continuum electrostatic
calculations. The atomic charges for the binuclear center,
successively with His291 protonated and deprotonated, cor-
respond to the Singh-Kollman charges calculated for all the
large Cu-Fe complexes. These ESP charges are given as
Supporting Information. The pKas were determined successively
including only segment A and within all four segments (A, B,
C, D) of the bacterial CcO. The protonation state of titratable
sites close to His291 could influence the acidity of imidazole.
Therefore, we calculated the pKa values successively with a
deprotonated and protonated Glu242. Located at the top of the
D-channel, Glu242 has been shown to be a key element in
proton translocation. This carboxylic group is thought to function
as a relay group shuttling protons from the D-channel to the
pocket of the binuclear center.6,60-66 In the crystal structure of
a bacterial CcO (i.e., PDB code 1M56, ref 6), the distance
between Glu242 carboxylic oxygen and His291 N2 is 11 Å
(see Figure 7). The protonation state of Glu242 throughout the
catalytic cycle has been the focus of many experimental studies
(see refs 67 and 68 and references therein). It has been shown
that this carboxylic group has an unusually high pKa of
approximately 9. However, it is not clear if it remains protonated
in all the significant catalytic steps.
The protonation state of the other titratable sites in the protein
corresponds to their standard state at pH ) 7. Subunits A and
D are neutral at pH ) 7. The few excess charges in subunits B
and C (and A when Glu242 is protonated) were neutralized by
protonation of solvent-exposed charged residues (e.g., Glu, Asp,
and Lys). These residues are located on the cytoplasmic side
of the membrane, and they are most likely neutralized by
counterions in vivo.
The atomic charges of the residues surrounding the binuclear
center correspond to the standard charges in the CHARMM22
parameter set.57 The Cu-Fe models used to derive the ESP-
fitted charges describe only a small part of the binuclear center.
Hence, all charges used in the continuum electrostatic calcula-
tions were rescaled to obtain the appropriate integer value of
the total charge. The dielectric constant inside the pocket
containing the ab initio-derived charges of the binuclear center
is set to QM ) 1. The protein reaction field was represented by
a dielectric constant p of 4, which is the value generally used
for protein calculations, and repeated with p ) 10. Recent
studies show that, because of water penetration, the polarizability
in the protein interior is described more accurately by a uniform
continuum dielectric with p between 8 and 12.69,70-72 Outside
of CcO, we set  ) 80 with an ionic strength of 150 mM.73
All geometry optimizations, energy, vibrational analysis, and
ESP charges calculations were performed with GAMESS-US.74
GB-CM3 calculations were performed with version 4.1 of
GAMESSPLUS.75 The MEAD package50 was used for the post-
SCF PB continuum electrostatic calculations.
3. Results
3.1. Choice of a Reference Acid-Base Couple. Continuum
solvation models commonly used to derive ¢Gsolv are able to
define very accurately the solvation energy of neutral species.
However, the electrostatic field of ionic species strongly perturbs
the solvent organization, particularly for solvent molecules in
the first solvation shell. These show very different properties
compared to the molecules in the bulk, hence can hardly be
described by a continuum solvation model. The discrepancy
between experiment and continuum electrostatic calculations of
ion solvation energies is expected to be in the range of 5
kcalâmol-1.30,43 To select an appropriate reference BH for the
calculation of the pKas within the relative scheme (Figure 5),
we tested the accuracy of the GB and PB solvation models in
the calculation of the solvation free energy for a few common
acid-base couples. The results are shown in Table 1. The largest
errors are associated with the solvation free energy of anions.
GB shows the largest deviation for the acetamide conjugate base.
The equilibrium MeNH3+ f MeNH2 shows the smallest
deviation for the total hydration free energy calculated with both
GB (j¢GsolGB - ¢Gsolexpj ) 0.606 kcalâmol-1) and PB (j¢GsolPB -
¢Gsol
exp
j ) 0.433 kcalâmol-1) solvation models. For this reason,
the couple MeNH3+/MeNH2 was chosen as reference. The
experimental pKa in bulk water for the deprotonation of
MeNH3+ is 10.66 (ref 40).
Figure 7. Binuclear center of cytochrome c oxidase (CcO) from PDB
structure 1M56. CuB is shown in teal, while Fe of hemea3 is shown in
orange. Hydrogen atoms were not included in the crystal structure.
TABLE 1: Solvation Free Energies in Bulk Water (E )
78.4) for Reference Equilibria (All Energies Are Expressed
in kcalÆmol-1)
acid/base ¢Gsol
GB(a) ¢Gsol
PB(b) ¢Gsol
exp(c) pKa
exp(c)
CH3CONH2/CH3CONH-1 -60.541 -69.098 -70.39 15.1
MeNH3+1/MeNH2 72.536 71.497 71.93 10.66
PhOH/PhO1- -58.333 -59.842 -64.68 9.99
pyridineH+1/pyridine 54.750 54.933 56.41 5.23
CH3COOH/CH3COO1- -67.240 -70.002 -70.60 4.76
PhNH3+1/ PhNH2 62.528 59.437 67.31 4.87
CH3CONH3+1/CH3CONH2 65.840 65.051 64.09 -0.66
mean abs. deviation maximum deviation
¢Gsol
GB ¢Gsol
PB Gsol
GB Gsol
PB
cations 2.809 2.048 6.452(d) 6.007(d)
anions 8.025 2.048 11.193(e) 4.335(g)
neutral 1.551 0.740 2.500(f) 1.866(h)
a Solvation free energy calculated with an SCF generalized Born
(GB) approach. b Solvation free energy calculated with a post-SCF
Poisson-Boltzmann (PB) approach. c Experimental values from ref 76.
d PhNH3+1. e CH3CONH-1. f MeNH2. g PhO1-. h PhNH2.
22634 J. Phys. Chem. B, Vol. 109, No. 47, 2005 Fadda et al.

3.2. pKa Values in Bulk Water. The relative scheme shown
in Figure 5 was tested successively on the deprotonation of
imidazole, which has an anionic conjugate base (i.e., C3N2H4
+ H2O f C3N2H3- + H3O+), and on the deprotonation of
imidazolium, which has a neutral conjugate base (i.e., C3N2H5+
+ H2O f C3N2H4 + H3O+). The experimental pKas reported
in the literature are 14 (refs 20-22) and 7.0 (ref 33) for
imidazole and imidazolium, respectively. The results are shown
in Table 2. For the deprotonation of imidazole, we obtain very
accurate results with GB (pKa ) 14.4), while the PB model
significantly overestimates the imidazole acidity, assigning a
pKa of 10.1. The opposite trend is shown in the deprotonation
of imidazolium, where GB slightly underestimates the pKa,
whereas PB shows very good accuracy.
The results obtained for imidazole as a ligand in different
Cu models are also shown in Table 2. For the deprotonation of
ImH in a Cu model with copper in its reduced state (i.e., CuI),
we obtain a pKa of 14.7 with the GB model, and 14.4 with the
PB model. The agreement between the two solvation methods
is most likely due to their similar level of accuracy in describing
the solvation of neutral species and cations (see Table 1). If we
compare these results with the ones obtained for the isolated
imidazole (pKaexp  14), it is clear that the presence of the
metal does not have any significant effect on the acidity.
Two Cu models with copper in its oxidized state (CuII) were
studied. In the first complex, CuII is bound to three imidazoles
rings and one hydroxyl, making the total charge of the acid
form equal to +1. In the other complex, the metal is ligated to
a water molecule as fourth ligand, so that the total charge of
the acid form is +2. The pKa of the imidazole ring in the
[ImH]3CuII-OH complex is very close to the one calculated
for the [ImH]3CuI-H2O despite the different oxidation state of
the metal. An increase in acidity is calculated for the [ImH]3CuII-
H2O complex. The results obtained with the GB model show a
decrease of 2.8 pKa units compared to the [ImH]3CuII-OH
complex. A less significant decrease of 0.9 pKa units is obtained
when using the PB solvation model. This change in acidity could
be due to the fact that the [ImH]3CuII-H2O complex is charged
in both protonated (+2) and deprotonated (+1) forms, whereas
both of the other two complexes have neutral conjugate bases.
The effect of the high-spin Fea3 group on the acidity of ImH
was investigated through the study of different Cu-Fe com-
plexes (see Figure 3). Each of the Cu-Fe models is designed
to reproduce a stage in the CcO catalytic cycle to assess the
possibility of a deprotonation of ImH linked to a proton pumping
event. Because complexes of CuI do not participate in pumping
events, in all of the Cu-Fe complexes studied, Cu is in its
oxidized state (i.e., CuII). The calculated pKa values are shown
in Table 3. The inclusion of the Fe-porphyrin results in an
increase of the imidazole pKa values, independently of the nature
of the intermetal ligands or the Fe redox state. The calculated
pKas are in the range of 14.8-16.6. The lowest pKa value
corresponds to the hypothetical ferromagnetic (S ) 3) model
of the O-state,
(the boxes around “Cu” and “Fe” indicate the presence of the
three imidazoles, and of NH3 and porphyrin, respectively), which
we are considering in our calculations in order to determine
the effect of magnetic coupling on the acidity of Cu-bound
imidazoles. The highest pKa value corresponds to the antifer-
romagnetic state (S ) 2).
The pKa values determined for the Cu- and Cu-Fe com-
plexes in bulk water can be used in eq 2 to obtain the proton
solvation free energy (¢GsolH
+). The results are shown in Table
4. All ¢Gsol data listed refer to the PB calculations for
consistency. The average of the ¢Gsol
H+
values derived from
TABLE 2: Calculated pKas in Bulk Water (E ) 78.4) of Isolated 5-Me-imidazole, Imidazolium, and of Cu-bound Imidazole
Ring in Various Minimal Cu Models (All Energies Are Expressed in kcalâmol-1)
acid ¢Hele ¢GTh ¢Gsol
GB(a) ¢Gsol
PB(b) pKa
GB pKa
PB
imidazole 359.175 -7.724 -58.039 -64.977 14.4 10.1(c)
imidazolium 234.416 -7.846 55.404 55.730 6.1 7.1(d)
CuI[(ImH)3 H2O]+1 281.271 -7.733 20.217 18.799 14.7 14.4
CuII[(ImH)3 OH]+1 276.688 -6.568 23.006 22.231 14.3 14.4
CuII[(ImH)3 H2O]+2 201.285 -7.179 95.225 96.972 11.5 13.5
Reference
MeNH3+1 224.082 -8.340 72.536 71.497 10.66(e)
a Solvation free energy calculated with an SCF generalized Born (GB) approach. b Solvation free energy calculated with a post-SCF Poisson-
Boltzmann (PB) approach. c The corresponding experimental pKa value is 14.0.20-22 d The corresponding experimental pKa value is 7.0.33 e Experimental
pKa from ref 76.
TABLE 3: Calculated pKa Values (E ) 78.4) of ImH in Various Cu-Fe models in Bulk Water (All Energies Are Expressed in
kcalâmol-1)
a Solvation free energy calculated with a post-SCF Poisson-Boltzmann (PB) approach. b The boxes around “Cu” and “Fe” stand to indicate the
presence of the three imidazoles and of NH3 and porphyrin, respectively.
Acidity of a Cu-Bound Histidine in Cytochrome c Oxidase J. Phys. Chem. B, Vol. 109, No. 47, 2005 22635

each equilibrium is -266.340 kcalâmol-1. The corresponding
standard deviation is 0.036 kcalâmol-1.
3.3. pKa Values in Cytochrome c Oxidase. As shown in
Figure 1, the CuB-hemea3 binuclear center is located in a pocket
deep inside subunit A of CcO. The effect of the proteic
environment on the acidity of the CuB-bound histidine is studied
by combining DFT with PB continuum electrostatics. Tables 5
and 6 list the results obtained via the thermodynamic cycle
in Figure 4, in combination with the ¢Gsol
H+ ) -266.340
kcalâmol-1 derived from the pKas in bulk water (see above).
The pKa values listed in Table 5 were obtained successively
from calculations including only subunit A and the full protein
(i.e., subunits A, B, C, and D in bacterial CcO), and with
deprotonated and protonated Glu242. Table 5 shows the results
obtained with p ) 4 to describe the protein reaction field, while
Table 6 shows the results obtained with p ) 10, a value which
takes into account the effect of water penetration in the protein
interior.69-71 Compared to the results in bulk water, the
consistent effect of the low dielectric proteic environment is to
shift the equilibrium toward higher pKa values. For antiferro-
magnetic coupled metals in the binuclear center, the pKa in CcO
ranges from 16.8 to 26.9 with p ) 4 and deprotonated Glu242.
Protonation of Glu242 results in a downshift of pKa values by
an average of 0.5 units when only subunit A is considered, and
by 2 units when the whole protein is included in the model
system (see Table 5). However, when the protein dielectric
constant is increased to p ) 10, the pKa values increase again
by 1.3 units (on average) for all complexes other than
whose pKa decreases by 1.7 units. This effect can be understood
in terms of the total charge of this complex, which, contrary to
the other states studied, has a neutral conjugate basis. Finally,
we note that pKa values of the same order of magnitude as the
ones calculated in the protein interior were obtained for the small
Cu models by continuum electrostatics in dimethyl sulfoxide
(DMSO,  ) 46.8). These data are included as Supporting
Information. This result is consistent with the previously
reported finding that a dielectric constant higher than 2-4 is
needed in continuum dielectric calculations in order to obtain
pKa values closer to the ones measured in proteins.77-79
4. Discussion
Within the framework of a recently proposed pumping
mechanism,14,15,19 the deprotonation of the CuB-bound H291 is
presumed to be triggered by the entry of a chemical proton in
the binuclear center, which converts a hydroxyl ligand into a
water molecule.14,15 According to this mechanism, H291 depro-
tonation takes place without dissociation of the N2-Cu bond,
which results in the formation of an imidazolate anion.
In the present work, various Cu- and Cu-Fe complexes were
modeled after the CcO binuclear center at different stages of
the catalytic cycle. The pKa values of the Cu-bound imidazole
in these complexes were determined successively in bulk water
and within CcO. The results obtained for the small Cu models
(Table 2) demonstrate that the effect of the oxidation state of
the Cu center does not influence the imidazole pKa. Indeed,
both [ImH]3CuI-H2O and [ImH]3CuII-OH, which bear the
same charge, but where the metal is in different redox states,
have very similar pKas. The latter are of the same order of
magnitude as the experimentally determined pKa of 14, corre-
sponding to the deprotonation of an isolated imidazole, and are
also in agreement with our value of 14.4 calculated with the
GB solvation model. The discrepancy of the corresponding pKa
determined with the PB approach is due to the difficulty intrinsic
to continuum solvation models in the description of anionic
species. Results from test calculations on common acid-base
TABLE 4: Proton Solvation Free Energy (¢GSolH
+) Derived
from Calculated pKa Values in Bulk Water (All Energies
Are Expressed in kcalâmol-1)
pKa(a) ¢Hgas(b) ¢Gsol(c) ¢Gsol
H+
7.1 220.250 55.730 -266.296
10.1 345.131 -64.977 -266.378
14.4 267.218 18.799 -266.375
14.4 263.800 22.231 -266.389
13.5 187.786 96.972 -266.344
15.0 220.260 66.535 -266.335
15.4 274.970 12.355 -266.320
14.8 216.544 69.955 -266.311
16.6 223.175 65.835 -266.367
15.8 219.906 67.932 -266.287
average -266.340
standard deviation 0.036
a pKa values (from ¢GsolPB) correspond to acid-base couples shown
in Tables 2 and 3. b ¢Hgas ) ¢Hele + ¢GTh + Htrans
H+ + ¢(pV) -
T[SH+]. c ¢Gsol ) ¢GsolA- - ¢GsolAH.
TABLE 5: Calculated pKas of the CuB-bound Imidazole
Ring of His291a
a Results are shown for deprotonated and protonated Glu242 suc-
cessively with only subunit A and with all four subunits (full protein),
with a protein dielectric constant  ) 4. All energies are expressed in
kcalâmol-1. b Binuclear center in CcO subunit A only, which is in turn
immersed in a continuum of  ) 80 with ionic strength 150 mM.
c Binuclear center surrounded by all four subunits (A, B, C, D) of CcO.
The protein is immersed in a continuum of  ) 80 with ionic strength
150 mM.
TABLE 6: Calculated pKas of the CuB-bound Imidazole
Ring of His291 with All Four Subunits (Full Protein), with a
Protein Dielectric Constant Ep ) 10a
a All energies are expressed in kcalâmol-1. b Binuclear center sur-
rounded by all four subunits (A, B, C, D) of CcO. The protein is
immersed in a continuum of  ) 80 with ionic strength 150 mM.
22636 J. Phys. Chem. B, Vol. 109, No. 47, 2005 Fadda et al.

couples, for which the solvation free energy is known experi-
mentally, show that the solvation of neutral species is very well
reproduced by both PB and GB methods (see Table 1). In this
particular case, the deviation in the solvation free energy
calculated for the imidazolate with PB is probably larger than
the one obtained for the GB method, and it is not counterbal-
anced by a similar deviation in the solvation free energy of the
acid. These results emphasize that, because pKa calculations are
extremely sensitive to small energy deviations, solvation free
energy contributions should be tested with different levels of
theory.
An increase in acidity is determined for the deprotonation of
[ImH]3CuII-H2O. This behavior can be interpreted in terms of
stabilization of charged species in bulk water. Here the acid
has a higher charge (+2) than that of the other complexes
studied. Its deprotonation leads to a charged conjugate base,
whereas both of the other two complexes have neutral conjugate
bases. One might consider [ImH]3CuII-H2O as produced by
the addition of a proton to the hydroxyl ligand of [ImH]3CuII-
OH. According to the proposed pumping mechanism, the
imidazole acidity should increase. Indeed, the pKa value in
[ImH]3CuII-H2O is 11.5 (according to the GB model), which
is significantly lower than that of [ImH]3CuII-OH (pKa ) 14.3),
whereas the corresponding pKa downshift calculated with the
PB solvation model is only of 0.9 units. Thus, despite the
relatively higher acidity of [ImH]3CuII-H2O, all Cu-imidazoles
complexes studied clearly favor the fully protonated form. This
is in agreement with the results of a computational study of
FeII-His pKa values in the Rieske protein active site.80 Results
show that the acidity of the His side chain bound to FeII is of
the same order of magnitude as that of free methyl-imidazole
(i.e., pKa values in the range of 11.3-12.8). Hence the effect
of the ligation to FeII on His acidity can be considered negligible,
and the complexes' fully protonated form is favored. Despite
the different stability of CuII- and FeII-His complexes,81 we
do also find that the pKa values of different CuII-His complexes
in bulk water are comparable to that of the free imidazole pKa
(i.e., approximately 14 pKa units).
Cu models are a very simplistic approximation of the CcO
binuclear center. We included the effect of the high-spin Fe-
hemea3 by studying larger Cu-Fe models with Fe in different
oxidation states. We can analyze the results within the frame-
work of the CcO catalytic cycle shown in Figure 2. The different
Cu-Fe models were designed to represent the various catalytic
steps connected to the deprotonation of the Cu-bound His291.
The first deprotonation is expected after the entry of a chemical
proton in F to form state Fp. The pKa in bulk water calculated
for a model of Fp, represented by
is 15.8, much higher than the pKa obtained for the corresponding
Cu model, [ImH]3CuII-H2O.
The next deprotonation is expected in connection with the
formation of the O-state, represented by
The corresponding pKa is 15.0. This indicates that the different
Fe redox states, consistently with the results obtained for the
various Cu redox states, do not cause a significant change on
the imidazole pKa. We also tried to study a different model of
the O-state, where the water molecule was ligated to Cu, while
the hydroxyl was ligated to Fe. However, the optimization of
the intermetal ligands leads inevitably to the complex
We studied two more states, both with Cu ligated to a
hydroxyl as a fourth ligand, and Fe with a free apical
coordination site. In the latter configuration, two oxidations
states of Fe were considered, Fe III and FeII. Because the
hydroxyl can function as a bridging ligand coordinating both
Cu and Fe, a ferromagnetic coupling between the two metals
can be envisaged, although it should be noted that only
antiferromagnetic states (S ) 2) were identified experimen-
tally.82 We modeled
with both ferromagnetic (S ) 3) and antiferromagnetic (S ) 2)
coupling to gauge the effect of spin coupling on the Cu-bound
imidazole pKa. The results in bulk water (Table 3) show that
the antiferromagnetic form
has a pKa slightly higher than the other complexes, 1.6 units
above the value determined for
which has the same total charge and oxidation state. The pKa
of the ferromagnetic form
is similar to that of
Given the limits of accuracy of the calculation, the magnetic
coupling does not seem to affect the imidazole acidity signifi-
cantly.
The pKa values obtained with the PB solvation model for
the small Cu models only differ by a few units (up to a
maximum deviation of 2.3 pKa units) from those obtained for
the larger Cu-Fe models in bulk water (see Tables 2 and 3).
While the Cu models were fully optimized in the protonated
and deprotonated form, for Cu-Fe models, only the protonated
form was relaxed within the hypothesis that deprotonation was
instantaneous. The agreement between the two sets of calcula-
tions indicates that structural modifications arising from geom-
etry optimization of the Cu-His complexes have a marginal
effect on the acidity of imidazole.
When we consider the Cu-Fe complexes as part of the active
site in the protein, the pKa of the CuB-bound imidazole ring of
H291 shifts to higher values in all relevant catalytic states (see
Tables 5 and 6). To gauge the effect of titration of protein side
chains on the acidity of H291, in this calculation, we have also
taken into account the protonation state of Glu242 (see Table
5). This group is the titratable group closest to His291, and it
is thought to play a role in the pumping mechanism.6,60-66 When
we consider only subunit A, the results obtained with a
protonated Glu242 are downshifted by only 0.6 pKa units (on
average) relative to the acidities determined when the Glu242
is deprotonated. When we consider all four subunits, this
Acidity of a Cu-Bound Histidine in Cytochrome c Oxidase J. Phys. Chem. B, Vol. 109, No. 47, 2005 22637

downshift is on average 2 pKa units. These pKa shifts correspond
to changes in the total solvation energy of only 0.8-2.8
kcalâmol-1, values which are within the limits of accuracy of
continuum solvation models.30,43 Thus the titration of Glu242
does not have a significant effect on the acidity of His291.
In two recent publications, Quenneville et al.36 and Popovic´
et al.37 studied the imidazole acidity in complexes similar to
ours. Results were reported for bulk water and in a uniform
dielectric of 4,36 as well as in the proteic environment through
continuum electrostatic calculations.37 Their results in bulk water
are similar to ours for [ImH]3CuIH2O. However, a pKa of 8.6
was calculated for the oxidized complex [ImH]3CuIIH2O, a much
lower value than the 11.5 (GB) and 13.3 (PB) we determined.
Because the protocol used to calculate the gas-phase energy
contribution is essentially the same, the discrepancy in this case
is due to relatively small differences (on the order of 2-5
kcalâmol-1) in the solvation energy contribution.
The most significant divergence pertains to the studies in low
dielectric and proteic environments. In a uniform continuum
dielectric of p ) 4, the authors report pKa values of -9.9 for
a CuII complex, and 9.8 for a CuI complex.36,37 The protein
reaction field was simulated also by including, within the
continuum dielectric of  ) 4, regions of  ) 80 corresponding
to protein cavities, which are presumably occupied by internal
water.37 By adding to this heterogeneous reaction field the
contribution of protein charges, Popovic´ et al.37 obtained pKa
values ranging between 2.1 and 4.0 for CuII complexes, and
17.4-19.7 for reduced CuI complexes. While water penetration
is likely to modulate the reaction field of the protein, the basis
for modeling internal hydration with dielectric boundaries in
continuum electrostatic models is unclear. First, choosing the
effective shape and the size of buried cavities is not a trivial
matter because the hydration state of cavities is modulated by
small changes in the structure and the polarity of the cavity
itself.83 Thus, it is likely that the hydration state of cavities
embedded in the interior of redox proteins such as CcO depends
on the redox state of the enzyme. In addition, evidence from
MD studies suggests that confinement significantly modifies
the relaxation dynamics of water molecules in pores and cavities
compared to that of bulk water.84-88 As a result, a dielectric of
80 is inadequate to represent confined water molecules. Finally,
the partition of the protein interior into many regions of
significantly different dielectric constants can lead to large errors
near dielectric boundaries. The reaction field contribution to
the free energy of ions in narrow pores was shown to be highly
sensitive to the proximity of dielectric boundaries between
regions of high and low dielectric.89 Recent studies suggest that
the effect of buried water on the protein reaction field may be
incorporated by applying a uniform dielectric constant between
8 and 12 throughout the protein interior.69-71 In the present
study, increasing the dielectric constant of the protein from p
) 4 to p ) 10 was found to cause the pKa values to increase
only by 1.4 units on average (see Tables 5 and 6). The only
exception is provided by
which, contrary to the other complexes, has a neutral conjugate
base (see Table 6).
Ultimately, the discrepancy between the results obtained by
Popovic´ et al.37 and by us can be related to three main factors.
First, in the previous studies, the pKa values were calculated
relative to the deprotonation of Me-imidazolium in bulk water
(pKa ) 6.6), even when the reaction takes place in an
environment different from bulk water, such as the uniform
dielectric with  ) 4, or the CcO protein interior. Yet, on the
basis of the thermodynamic cycle depicted in Figure 5, it is
clear that, in order to eliminate the dependence from the proton
solvation free energy, both equilibria have to take place in the
same solvent. Because the expected pKa value of imidazolium
in a low dielectric is presumably higher than 6.6 (e.g., a pKa
value of 10.2 was determined for His63 in thiolsubtilisin
Carlsberg90), the data reported by Popovic´ et al.37 are likely to
be artificially downshifted. Furthermore, as indicated above,
dielectric boundaries separating regions of low dielectric ( )
4 and 1) from regions of high dielectric ( ) 80) may cause
large fluctuations in the reaction field, which could be critical
to the accuracy of acidity calculations if such boundaries are
located close to the deprotonation site in the binuclear center.
Finally, another significant difference between the present study
and that of Popovic´ et al.37 is that the ESP charges used in their
PB calculation for the whole protein were calculated only on
the small CuB complex, which does not include hemea3. As a
result, their electrostatic potential may be inappropriate to
describe the binuclear center, especially in the intermetal region.
The results that we obtained have important implications for
the proposed pumping mechanism. Within the hypothesis of
the CuB-bound H291 imidazole ring as the loading site of the
CcO proton pump, we would expect a strong thermodynamic
coupling between the redox centers (CuB and/or hemea3) and
the protonation site. This entails a significant dependence of
the pKa of the proton loading site on the redox state of the
binuclear center. The calculations performed on Cu models show
no significant dependence of pKa on the Cu oxidation states.
When the pumping events take place along the catalytic cycle
(see Figure 2), Cu is always present as CuII. The tests conducted
in bulk water on larger Cu-Fe complexes, which reproduce
more closely the binuclear center, do not support deprotonation
of the Cu-bound imidazole. The lower dielectric of the proteic
environment has the effect of shifting the already high pKas to
even higher values, indicating a considerable destabilization of
the CuB-bound imidazolate produced by deprotonation within
the protein. It is likely that this destabilization would be less
dramatic if we had included explicit water molecules in the
hydrophobic pocket of the active site. Nevertheless, this
refinement of the model would not change our conclusion
regarding the absence of thermodynamic coupling between the
redox states of the binuclear center and the pKa of the CuB-
bound imidazole.
It should be noted that a study of pKa, an equilibrium quantity,
falls short of the level of detail required to encompass the proton
pumping mechanism because pumps work out of equilibrium.
The inclusion of transient dynamic events that are likely to play
a role in the proton pumping mechanism is beyond the scope
of this study. With this caveat in mind, the significance of the
pKa values obtained in the present study resides in the
considerable energetic cost associated with deprotonation of the
Cu-bound imidazole. For that reason, our findings are incon-
sistent with a mechanism in which H291 is the key pumping
element of the enzyme.
5. Conclusions
In this work, we calculated pKa values for Cu-bound
imidazoles in various Cu- and Cu-Fe complexes, with the
purpose of assessing the validity of a recently proposed proton
pumping mechanism for CcO, based on the deprotonation of a
Cu-bound histidine side chain.14,15 Our approach is based on
the study of increasingly more detailed models of CcO binuclear
22638 J. Phys. Chem. B, Vol. 109, No. 47, 2005 Fadda et al.

center, which allows us to analyze successively the effect of
the redox states of Cu, of Cu and Fe, and of the environment
surrounding the binuclear center (i.e., bulk water vs enzyme
interior). According to the proposed catalytic scheme,14,15 His291
undergoes deprotonation whenever a chemical proton enters the
binuclear center and converts the hydroxyl ligated to one of
the metals (Cu or Fe) into water. For this mechanism to work,
the imidazole ring pKa must change considerably depending on
the different catalytic steps. Hence, a thermodynamic coupling
between the binuclear center redox states and the CuB-bound
imidazole ring acidity is to be expected.
The experimental pKa value for the deprotonation of imida-
zole, which produces an anionic imidazolate, is approximately
14.20-22 Our calculations in bulk water on small Cu complexes
do not show any significant change of the pKa due to the ligation
of the imidazole to Cu. Furthermore, the results obtained for
CuI and CuII complexes bearing the same total charge do not
show any significant dependence on the redox state of the metal.
The study of larger complexes, which include an Fe-
porphyrin with Fe in different redox states (i.e., II, III, and IV),
representing the hemea3 at various stages of the catalysis, reveals
a small increase in the pKa of imidazole. Still, no significant
dependence of the imidazole pKa on the oxidation state of the
binuclear center was determined.
Finally, we considered the effect of the protein environment
by combining DFT calculations with PB continuum electrostat-
ics. The effect of the proteic environment was to increase further
the pKa values in all complexes, independently of the redox
state of the metals. The calculations in protein required the use
of an absolute scheme, which takes into account explicitly the
contribution proton solvation free energy (¢GsolH+). The experi-
mental values attribued to ¢Gsol
H+
vary by as much as 12
kcalâmol-1. Therefore, for consistency with the calculations in
bulk water, we derived our estimate of ¢Gsol
H+
. We obtained
¢Gsol
H+ ) -266.340 kcalâmol-1, in excellent agreement with
accurate measurements reported in the literature.31 This result
is significant, especially in view of the difficulties intrinsic in
the experimental and computational assessment of ¢Gsol
H+
and
of the consequent uncertainty on the ¢Gsol
H+
value. Our estimate
supports the finding that ¢Gsol
H+ is much more negative than
originally thought.31
All the pKas obtained in bulk water as well as in protein are
considerable. The values range between 14.8 and 16.6 in bulk
water, and between 14.9 and 26.9 in the protein. Continuum
calculations based on static molecular models inherently neglect
transient events that are required for kinetic control of the
process, giving rise to the directional translocation of protons
against an electromotive force. Nevertheless, the present results
indicate that the energy needed to deprotonate the Cu-bound
histidine is very high, which is inconsistent with its previously
proposed role as a key proton loading site for the pump
mechanism.
The present study should be seen as a necessary first step in
characterizing important aspects of the chemistry and equilib-
rium properties of the enzyme. Future efforts will focus on
examining the interplay between equilibrium properties, such
as redox and hydration states, and electron and proton transport
events.
Acknowledgment. We are grateful to Dr. Ching-Hsing Yu
for his help with the computational resources and to the Center
for Computational Biology for a generous allocation of computer
time. We thank Dr. G. C. Shields and Dr. V. Guallar for helpful
discusions and Dr. D. Bashford for providing us with the
program MEAD. The Hospital for Sick Children Training Center
and the Canadian Institutes for Health Research (CIHR operating
grant M0P43949) are gratefully acknowledged for financial
support. R.P. is a CRCP chairholder.
Supporting Information Available: Cartesian coordinates
and electrostatic potential fitted charges for all complexes. pKa
values in DMSO for Cu models. This material is available free
of charge via the Internet at http://pubs.acs.org.
References and Notes
(1) Wikstro¨m, M. Nature 1977, 273, 271.
(2) Gennis, R. B. Front. Biosci. 2004, 1, 581.
(3) Saraste, M. Science 1999, 283, 1488.
(4) Ferguson-Miller, S.; Babcock, G. T. Chem. ReV. 1996, 96, 2889.
(5) Wikstro¨m, M. Curr. Opin. Struct. Biol. 1998, 8, 480.
(6) Svensson-Ek, M.; Abramson, J.; Larsson, G.; To¨rnroth, S.; Brzez-
inski, P.; Iwata, S. J. Mol. Biol. 2002, 321, 329.
(7) Iwata, S.; Ostermeier, C.; Ludwig, B.; Michel, H. Nature 1995,
376, 660.
(8) Tsukihara, T.; Aoyama, H.; Yamashita, E.; Tomizaki, T.; Yamagu-
chi, H.; Shinzawa-Itoh, K.; Nakashima, R.; Yaono, R.; Yoshikawa, S.
Science 1995, 269, 1069.
(9) Gennis, R. B. Biochim. Biophys. Acta 1998, 1365, 241.
(10) Wikstro¨m, M.; Jasaitis, A.; Backgren, C.; Puustinen, A.; Verkhovsky,
M. I. Biochim. Biophys. Acta 2000, 1459, 514.
(11) Wikstro¨m, M. Biochim. Biophys. Acta 2000, 1458, 188.
(12) Wikstro¨m, M.; Verkhovsky, M. I. Biochim. Biophys. Acta 2002,
1555, 128.
(13) Wikstro¨m, M.; Verkhovsky, M. I.; Hummer, G. Biochim. Biophys.
Acta 2003, 1604, 61.
(14) Popovic´, D. M.; Stuchebrukhov, A. A. J. Am. Chem. Soc. 2004,
126, 1858.
(15) Popovic´, D. M.; Stuchebrukhov, A. A. FEBS Lett. 2004, 566, 126.
(16) Gennis, R. B. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 12747.
(17) Michel, H. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 12819.
(18) Olkhova, E.; Hutter, M. C.; Lill, M. A.; Helms, V.; Michel, H.
Biophys. J. 2004, 86, 1873.
(19) Das, T. K.; Gomes, C. M.; Teixeira, M.; Rousseau, D. L. Proc.
Natl. Acad. Sci. U.S.A. 1999, 96, 9591.
(20) Walba, H.; Isensee, R. W. J. Am. Chem. Soc. 1955, 77, 5488.
(21) Yegil, G. Tetrahedron 1967, 23, 2855.
(22) Cleland, W. W.; Frey, P. A.; Gerlt, J. A. J. Biol. Chem. 1998, 273,
25529.
(23) Iwaki, M.; Yakovlev, G.; Hirst, J.; Osyczka, A.; Dutton, P. L.;
Marshall, D.; Rich, P. P. Biochemistry 2005, 44, 4230.
(24) Lodi, P. J.; Knowles, J. R. Biochemistry 1991, 30, 6948.
(25) Yoshikawa, S.; Shinzawa-Itoh, Nakashima, R.; Yaono, R.; Ya-
mashita, E.; Inoue, N.; Yao, M.; Fei, M. J.; Peters Libeu, C.; Mizushima,
T.; Yamaguchi, H.; Tomizaki, T.; Tsukihara, T. Science 1998, 280, 1723.
(26) Proshlyakov, D. A.; Pressler, M. A.; Babcock, G. T. Proc. Natl.
Acad. Sci. U.S.A. 1998, 95, 8020.
(27) Proshlyakov, D. A.; Pressler, M. A.; DeMaso, C.; Leykam, J. F.;
DeWitt, D. L.; Babcock, G. T. Science 2000, 290, 1588.
(28) Iwaki, M.; Puustinen, A.; Wikstro¨m. M.; Rich, P. R. Biochemistry
2004, 43, 14370.
(29) Koch, W.; Holthausen, M. C. A Chemist’s Guide to Density
Functional Theory, Wiley-VCH: Weinheim, 2000.
(30) Chen, J. L.; Noodleman, L.; Case, D. A.; Bashford, D. J. Phys.
Chem. 1994, 98, 11059.
(31) Liptak, M. D.; Shields, G. C. J. Am. Chem. Soc. 2001, 123, 7314.
(32) Reiss, H.; Heller, A. J. Phys. Chem. 1985, 89, 4207.
(33) Lim, C.; Bashford, D.; Karplus, M. J. Phys. Chem. 1991, 95, 5610.
(34) Marcus, Y. Ion SolVation, Wiley: New York, 1985.
(35) Tawa, G. J.; Topol, I. A.; Burt, S. K.; Caldwell, R. A.; Rashin, A.
A. J. Chem. Phys. 1998, 109, 4852.
(36) Quenneville, J.; Popovic´, D. M.; Stuchebrukhov, A. A. J. Phys.
Chem. B 2004, 108, 18383.
(37) Popovic´, D. M.; Quenneville, J.; Stuchebrukhov, A. A. J. Phys.
Chem. B 2005, 109, 3616.
(38) Atkins, P. Physical Chemistry, 6th ed.; W. H. Freeman & Co.: New
York, 1997.
(39) Chipman, D. M. J. Phys. Chem. A 2002, 106, 7413.
(40) Pliego, J. R., Jr.; Riveros, J. M. Chem. Phys. Lett. 2000, 332, 597.
(41) Brodskaya, E.; Lyubartsev, A. P.; Laaksonen, A. J. Phys. Chem. B
2002, 106, 6479.
(42) Toth, A. M.; Liptak, M. D.; Phillips, D. L.; Shields, G. C. J. Chem.
Phys. 2001, 114, 4595.
Acidity of a Cu-Bound Histidine in Cytochrome c Oxidase J. Phys. Chem. B, Vol. 109, No. 47, 2005 22639

(43) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. 1996,
100, 19824.
(44) Stevens, P. J.; Devlin, J. F.; Chabalowski, C. F.; Frisch, M. J. J.
Phys. Chem. 1994, 98, 11623.
(45) Stevens, W. J.; Basch, H.; Krauss, M.; Jasien, P. Can. J. Chem.
1992, 70, 612.
(46) Cundari, T. R.; Stevens, W. J. J. Chem. Phys. 1993, 98, 5555.
(47) Hay, P. J.; Wadt, W. R. J. Chem. Phys. 1984, 82, 299.
(48) Winget, P.; Thompson, J. D.; Xidos, J. D.; Cramer, C. J.; Truhlar,
D. G. J. Phys. Chem. A 2002, 106, 10707.
(49) Thompson, J. D.; Cramer, C. J.; Truhlar, D. G. J. Comput. Chem.
2003, 24, 1291.
(50) Bashford, D. Lect. Notes Comput. Sci. 1997, 1343, 233; Scientific
Computing in Object-Oriented Parallel EnVironments; Ishikawa, Y.,
Oldehoeft, R. R., Reynders, J. V. W., Tholburn, M. Eds.; ISCOPE97;
Springer: Berlin, 1997.
(51) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1984, 5, 129.
(52) Sigfridsson, E.; Ryde, U. J. Comput. Chem. 1998, 19, 377.
(53) Bayly, C. I.; Cieplack, P.; Cornell, W. D.; Kollman, P. A. J. Phys.
Chem. 1993, 97, 10269.
(54) Breneman, C. M.; Wiberg K. B. J. Comput. Chem. 1990, 11, 361.
(55) Chirlian, L. E.; Francl, M. M. J. Comput. Chem. 1987, 8, 894.
(56) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110,
165.
(57) MacKerell, J. A. D.; Bashford, D.; Bellot, M.; Dunbrack, J. R. L.;
Evanseck, J. D.; Field, M. J.; Fisher, G.; Gao, J.; Guo, H.; Ha, S.; Joseph-
McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.;
Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodholm, B.; Reiher, I. W. E.;
Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe,
M.; Wio´rkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. 1998,
102, 3586.
(58) Bondi, A. J. Phys. Chem. 1964, 68, 441.
(59) Baik, M.; Friesner, R. A. J. Phys. Chem. A 2002, 106, 7407.
(60) Pome`s, R.; Hummer, G.; Wikstro¨m, M. Biochim. Biophys. Acta:
Bioenerg. 1998, 1365, 255.
(61) Hellwig, P.; Rost, B.; Kaiser, U.; Ostermeier, C.; Michel, H.;
Mantele, W. FEBS Lett. 1996, 385, 53.
(62) Puustinen, A.; Bailey, J. A.; Dyer, R. B.; Mecklenburg, S. L.;
Wikstro¨m, M.; Woodruff, W. H. Biochemistry 1997, 36, 13195.
(63) Adelroth, P.; Svensson-Ek, M.; Mitchell, D. M., Gennis, R. B.;
Brzezinski, P. Biochemistry 1997, 36, 13824.
(64) Junemann, S.; Meunier, B.; Fisher, N.; Rich, P. R. Biochemistry
1999, 38, 5248.
(65) Rich, P. R.; Breton, J.; Junemann, S.; Heathcote, P. Biochim.
Biophys. Acta: Bioenerg. 2000, 1459, 475.
(66) Adelroth, P.; Karpefors, M.; Gilderson, G.; Tomson, F. L.; Gennis,
R. B.; Brzezinski, P. Biochim. Biophys. Acta 2000, 1459, 533.
(67) McMahon, B. J.; Fabian, M.; Tomson F.; Causgrove, T.; Bailey,
J. A.; Rein, F. N.; Dyer, R. B.; Palmer, G.; Gennis, R. B.; Woodruff, W.
H. Biochim. Biophys. Acta 2004, 1655, 321.
(68) Gennis, R. B. FEBS Lett. 2003, 555, 2.
(69) Damjanovic´, A.; Garcı´a-Moreno, B.; Lattman, E. E.; Garcı´a, A. E.
Proteins 2005, 60, 433.
(70) Archontis, G.; Simonson, T. Biophys. J. 2005, 88, 3888.
(71) Dwyer, J. J.; Gittis, A. G.; Karp, D. A.; Lattman, E. E.; Spencer,
D. S.; Stites, W. E.; Garcı´a-Moreno, B. Biophys. J. 2000, 79, 1610.
(72) Thompson, J. D.; Cramer, C. J., Truhlar, D. G. J. Phys. Chem. A
2004, 108, 6532.
(73) Bashford, D.; Karplus, M. Biochemistry 1990, 29, 10219.
(74) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.;
Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.;
Su, S. J.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem.
1993, 14, 1347.
(75) Pu, J.; Thompson, J. D.; Xidos, J. D.; Li. J.; Zhu, T.; Hawkins, G.
D.; Chuang, Y.-Y.; Fast, P. L.; Liotard, D. A.; Rinaldi, D.; Gao, J.; Cramer,
C. J.; Truhlar D. G. GAMESSPLUS, version 4.1; University of Minnesota:
Minneapolis, 2003.
(76) Pliego, J. R., Jr.; Riveros, J. M. Phys. Chem. Chem. Phys. 2002, 4,
1622.
(77) Dwyer, J. J.; Gittis, A. G.; Karp, D. A.; Lattman E. E.; Spencer,
D. S.; Stites, W. E.; Garcı´a-Moreno E. B. Biophys. J. 2000, 79, 1610.
(78) Antosiewicz, J.; McCammon, J. A.; Gilson, M. K. J. Mol. Biol.
1994, 238, 415.
(79) Antosiewicz, J.; McCammon, J. A.; Gilson, M. K. Biochemistry
1996, 35, 7819.
(80) Ullmann, G. M.; Noodleman, L.; Case, D. A. J. Biol. Inorg. Chem.
2002, 7, 632.
(81) Sjo¨berg, S. Pure Appl. Chem. 1997, 69, 1549.
(82) Ghilardi, R. A.; Huang, H. W.; Moe¨nne-Laccoz, P.; Stasser, J.;
Balckburn, N. J.; Woods, A. S.; Cotter, R. J.; Incarvito, C. D.; Rheingold,
A. L.; Karlin, K. D. J. Biol. Inorg. Chem. 2005, 10, 63.
(83) Vaitheeswaran, S.; Yin, H.; Rasaiah, J. C.; Hummer, G. Proc. Natl.
Acad. Sci. U.S.A. 2004, 101, 17002.
(84) Allen, T. W.; Kuyucak, S.; Chung, S. H. J. Chem. Phys. 1999,
111, 7985.
(85) Tieleman, D. P.; Berendsen, H. J. C. Biophys. J. 1998, 74, 2786.
(86) Green, M. E.; Lu, J. J. Phys. Chem. B 1997, 101, 6512.
(87) Sansom, M. S. P.; Smith, G. R.; Adcock, C.; Biggin, P. C. Biophys.
J. 1997, 73, 2404.
(88) Pome`s, R.; Roux, B. Biophys. J. 2002, 82, 2304.
(89) Chakrabarti, N.; Roux, B.; Pome`s, R. J. Mol. Biol. 2004, 343, 493.
(90) Kahyaoglu, A.; Jordan, F. Protein Sci. 2002, 11, 965.
22640 J. Phys. Chem. B, Vol. 109, No. 47, 2005 Fadda et al.