12162
J. Phys. Chem. B 2006, 110, 12162-12166
Improved Density Functional Theory/Electrostatic Calculation of the His291 Protonation State in Cytochrome c Oxidase: Self-Consistent Charges for Solvation Energy Calculation D. V. Makhov, Dragan M. Popovic´ , and Alexei A. Stuchebrukhov* Department of Chemistry, UniVersity of California, DaVis, One Shields AVenue, DaVis, California 95616 ReceiVed: February 9, 2006; In Final Form: April 19, 2006
The protonation state of His291 in cytochrome c oxidase (CcO), a ligand to the CuB center of the enzyme, has been recently studied in this group by using combined density functional theory (DFT)/electrostatic (QM/ MM) calculations. On the basis of these calculations, a model of the proton pumping mechanism of CcO has been proposed. Due to certain technical difficulties, the procedure used in the previous calculation to find partial atomic charges of the QM system for the solvation energy evaluation was not entirely satisfactory; i.e., it was not self-consistent. Here, we describe a procedure that resolves the problem and report on the improved calculations of the protonation state of the His residue. The new procedure fits the protein and reaction field potentials in the region of the QM system with artificial point charges placed on a surface of a sphere surrounding the QM system and a few charges inside the sphere and allows one to perform DFT calculations that involve an inhomogeneous dielectric environment in a self-consistent way. The procedure improves the accuracy of calculations in comparison with previous work. The improved results show, however, that although the absolute energies change significantly the relative energies of the protonated and deprotonated states of His291 remain close to the previously reported ones and therefore do not change significantly the pKa values reported earlier. Therefore, our new improved calculations support for the proposed His291 model of the CcO pump.
1. Introduction Cytochrome c oxidase (CcO)1-7 is a terminal enzyme in the respiratory electron transport chain of eukaryotic organisms and some aerobic bacteria. CcO catalyzes the reduction of oxygen to water and uses the energy of this reaction to pump protons across the inner mitochondrial membrane, a process that creates the membrane electrochemical proton gradient that drives the synthesis of ATP.8-12 The structure of CcO has been solved in the past decade for various organisms;13-17 however, the exact molecular mechanism of proton pumping is still mostly unknown.1-7,18,19 To identify the mechanism of proton pumping by CcO, a number of density functional theory (DFT)/electrostatic calculations have been recently performed in our group20-23 on the redox dependence of the protonation state (i.e., calculation of the pKa values) of CcO. These calculations point to a possible key role played by His291, one of the ligands to the CuB center of the enzyme, which seems to change protonation state depending on the redox state of the CuB center. These results suggest a possible mechanism of proton pumping in CcO.22 In combined ab initio/electrostatic methods, ab initio quantum chemical calculations are performed for the part of the protein (QM system) that actively participates in the reaction and includes the catalytic Fea3-CuB binuclear center and their ligands (Figure 1), while the rest of protein is treated as a polarizable dielectric environment with point charges. The pKa evaluation involves both quantum energies of the QM system and the solvation energy of the QM system in the protein environment. This dielectric environment is highly nonuniform * Author to whom correspondence should be addressed. Fax: 530-7528995. E-mail:
[email protected].
Figure 1. Surface rendering of the QM system representing the binuclear center of bovine CcO.16 Imidazole-Fea3(IV)-oxo porphyrin and the CuB center consisting of the three methylimidazoles and a ligated H2O molecule are displayed inside of the QM cavity. The two metal atoms, Fea3 and CuB, are shown as large spheres; the protonatable site of His291 is shown with a smaller sphere.
due to the presence of water cavities inside of the protein. These cavities have been shown to significantly affect the pKa value of His291.20,23 Unfortunately, the inclusion of water cavities into calculations presents significant difficulties. Not only the proper dielectric
10.1021/jp0608630 CCC: $33.50 © 2006 American Chemical Society Published on Web 05/27/2006
Improved DFT/Electrostatic Calculations for CcO properties of the cavities are unknown, but also the standard quantum chemistry codes typically would not provide the possibility to perform calculations in a nonuniform dielectric environment. This has been the case in our previous calculations20,23 that utilized the Jaguar program.24 Due to this technical difficulty, in our previous work, the DFT calculations for the QM system were performed in a uniform dielectric environment with ) 4 (which corresponds to a “dry” protein). The electrostatic potential (ESP) fitted atomic point charges for the QM system found in these calculations were then passed to a Poisson-Boltzmann (PB) solver, and the reaction field caused by these charges in a nonuniform environment was calculated. The interaction of the ESP atomic charges of the QM system with the polarization reaction field of the protein determines the solvation energy of the QM system in the protein environment. Thus the procedure used was not self-consistent. By using this procedure, we neglected both the effect of dielectric inhomogeneity of the environment and the protein charges on the distribution of the electron density in the QM system and corresponding ESP charges. The assumption was that the effect is likely to be equal in the protonated and deprotonated states of His291 and, therefore, does not affect the pKa values, which depend on the difference between solvation energies. Nevertheless, the question remained open about the legitimacy of this approach. To address this problem, in this paper we perform DFT/ electrostatic calculations for CcO using a new approach. We have developed a procedure to fit the electrostatic potential given by the PB solver to a field from the set of point charges. Because most standard quantum chemistry packages, including Jaguar, allow the inclusion of additional point charges, this makes it possible to perform ab initio calculations for the QM system in complex dielectric environments in a simple computational procedure. The plan of the rest of the paper is as follows. The details of the new computational procedure are described in section 2, while section 3 presents the results of the improved calculation for CcO and discusses the comparison of the new results with the previously published data.20 Concluding remarks are given in the last section. 2. Method of Calculations The QM system occupies a specific volume, a cavity inside of the protein. Inside of this volume, we set the dielectric constant to ) 1. The potential V acting on the QM system due to the surrounding protein medium can be divided into two parts. The first part is called reaction field; it represents a polarization response of the surrounding dielectric medium to the electrostatic field (charges) of the QM system. The second part is the field of the atomic protein charges, the protein field. As the reaction field, the protein field is nontrivial, because the Coulomb interaction of the protein charges and charges of the QM system is modified by the dielectric response of the medium. Thus both parts involve polarization response of the protein medium. Both potentials are found by solving the PB equation in the proper dielectric environment. However, to perform ab initio calculations in this environment, one should find a way to pass the potential given by the PB solver to the quantum chemistry code, because the distribution of electron density of the QM system, and its corresponding atomic charges, will depend on the “external” protein field acting on the QM system. The straightforward way to do this is to add this potential on the grid used by a DFT code when calculating the exchange-correlation potential.25,26 However, this requires a
J. Phys. Chem. B, Vol. 110, No. 24, 2006 12163 serious modification of the original DFT code, which is not always possible for standard ab initio packages (at the very least, one needs access to the original code, which is not always readily available). We propose an alternative way where the reaction and protein field potentials are fitted by the sum of point charge potentials. Practically, this is very convenient, because most of the standard quantum chemistry packages allow the addition of point charges to the system in a vacuum; thus this procedure would allow us to perform DFT calculations in the external potential given by the PB solver. Once this is implemented, a self-consistent procedure for doing such calculations can be readily realized. The reaction and protein field potentials are created by the permanent or induced charges outside of the cavity occupied by the QM systems. Thus, inside the QM cavity, the potential V must satisfy the Laplace equation
∆V ) 0
(1)
which has the following solutions in spherical coordinates
V(r,θ,φ) ) rlYlm(θ,φ)
(2)
where l g 0 and -l e m e l (we do not consider here the second set of solutions that is singular at the origin). Functions Ylm, known as spherical harmonics, are defined as
Ylm(θ,φ) )
{
Pml (cos θ) cos mφ 0 e m e l Pml (cos θ) sin mφ -l e m < 0
(3)
where Pml are Legendre functions. Inside the cavity, these solutions form a complete basic set consisting of the functions l V sph k ) r Ylm(θ,φ)
(4)
where k ) {l,m} is introduced to shorten notation. Thus, the potential acting on the QM system then can be represented as a sum
V)
sph ∑k asph k Vk
(5)
with coefficients ak to be determined from the fitting procedure. It can be easily shown27 that the basic functions (eq 4) represent an electrostatic potential inside a charged sphere of radius R with the charge distribution
Pk(θ,φ) )
2l + 1 l-1 R Ylm(θ,φ) 4π
(6)
Thus, if we consider a sphere enveloping our QM cavity (we call it the QM sphere) with the charge distribution
F)
∑k asph k Fk
(7)
the potential inside the cavity caused by this charged sphere will be identical to that given by the PB solver. In the case of the reaction field, the potential is usually smooth, and the small-length-scale components of the field can be neglected. This means that the high-frequency terms (large l) of the summation in eq 5 can usually be cut off at some reasonable l determined by the required accuracy of the calculations. Then, if we choose a spherical grid where the distances between grid points are smaller than the cutoff length scale, then we can replace the smooth charge density with a set of point charges on this grid. Then, the potential due to these
12164 J. Phys. Chem. B, Vol. 110, No. 24, 2006
Makhov et al.
point charges will reproduce with reasonable accuracy the reaction field of the nonuniform dielectric environment. Fitting Procedure. Reaction Field. Our fitting procedure for the reaction field therefore consists of the following steps. First of all, we choose a QM sphere enveloping the QM cavity (Figure of Nsph points on this sphere 1) and set up a uniform grid rsph j using a procedure described in ref 28. The radius of the sphere should not be too large to avoid extremely high charges for high-order harmonics; on the other hand, one should provide sufficient distance between the sphere and the QM system, at least larger than the distance between grid points on the sphere. Because the grid is uniform, the charge density Fk would correspond to the following set of point charges on the grid
Qsph kj )
4πR2 Fk(rsph j ) Nsph
(8)
Then, we set up a second grid Ri inside the QM cavity on which the potential will be fitted. To minimize fitting errors for both the potential and the electric field on atoms of the QM system, we use a nonregular grid that includes all atomic positions and sets of six additional points around each atom center shifted by 0.25 Å in each direction. We then calculate the potential Vsph ik generated on grid points Ri by the point charges on the QM sphere defined in eq 8, which will be used as a basis in potential fitting. By design, this potential is close to the spherical harmonics in eq 4; nevertheless, for better precision it should be calculated directly as
V sph ik
)
∑j
Qsph kj |Ri - rsph j |
(9)
To fit the reaction field potential, Vrf we now solve in the leastsquares sense the following overdetermined system of linear equations
Vrf(Ri) )
˜ sph ∑k Vsph ik a k
(10)
)
∑k
Qsph kj
a˜ sph k
V(Ri) )
∑j Aijqatj
(12)
where
1 |Ri - ratj |
Aatij )
(13)
In the SVD procedure, the (M × N) matrix ||Aatij || is presented as
Aatij )
∑k VatikwkQatkj
(14)
where wk’s are singular values of ||Aatij || and ||Vatik|| and ||Qatkj|| are orthogonal matrices
∑i Vatik Vatil ) ∑j Qatkj Qatlj ) δkl
(15)
The (M × M) matrix ||Vatik|| is square; its columns form a complete ortho-normal basis set of vectors for the expansion of potential V on M grid points in the QM system. The (M × N) matrix ||Qatkj|| is such that its M rows form a set of M orthonormal basis vectors for representing the distribution of charges. Thus, the potential can be expanded as M
which gives us contributions a˜ k from different basis functions. The tilde in the notation is used here to reflect the fact that the basis functions are slightly different from the spherical harmonics. Then, the point charges on the QM sphere can be found by linear transformation
qsph j
QM sphere. (These additional charges are typically required only for the protein field.) The procedure of fitting the protein field then consists of two stepssfinding appropriate values of the point charges, to reproduce the most singular part of the potential, and then using spherical harmonics to fit the remaining, nonsingular (i.e., smooth) part of the protein field. The point charges are chosen by utilizing the singular value decomposition (SVD) procedure29 as follows. The electrostatic potential created by N charges qatj on M grid points Ri in the QM region is given by the equation
(11)
Protein Field. A slightly different procedure is required for fitting the protein field. Some of the atoms of the protein are very close to the QM system (e.g., some of these atoms are bonded to the QM atoms). Such atoms can create a significant potential on a few (practically, one or two) atoms of the QM system, much larger than on other atoms. In other words, the protein field in the QM region is not a smooth function. An accurate fitting of such a potential by spherical harmonics is impossible unless the high-frequency harmonics (i.e., large l) are used. This is not convenient, because it would require too many charge grid points on the QM sphere; moreover, charge density increases with l and can become extremely large for high-frequency harmonics causing severe computational errors. This problem can be avoided by adding to spherical harmonics additional point charges qatj placed on the positions ratj of protein atoms that are external to the QM cavity and inside the
∑ aatk Vatik k)1
V(Ri) )
(16)
and the charges as M
qatj )
batk Qatkj ∑ k)1
(17)
(Notice that the basis set for N charges (N g M) is not complete but, however, is sufficient to represent M values of the potential.) Let the potential V(Ri) be given on the grid points of the QM system; its expansion coefficients are found from eq 16 using the orthogonality condition of the matrix ||Vatik|| M
aatk )
V(Ri)Vatik ∑ i)1
(18)
As it follows from the form of the matrix ||Aatij || eq 14, one possible set of charges that satisfies eq 12 is given by eq 17, with coefficients batk given by
batk
aatk ) wk
(19)
Improved DFT/Electrostatic Calculations for CcO
J. Phys. Chem. B, Vol. 110, No. 24, 2006 12165
Figure 2. Singular values of the matrix ||Aatij || (see text) for the case of the protonated oxidized state of the binuclear center.
Formally, the above relations, eqs 16-19, give an exact solution to eq 12; however, there is a problem. As it appears in eq 19, for small or zero singular values wk, the charges qatj , given by eq 17, become unrealistically large. And even if they faithfully reproduce the potential at the grid points, between the grid points the potential generated by such charges is likely to be unrealistically large too, because the large charges of both signs while compensating for the potential at the grid points will not do so between the grid points. This underscores the general problem that it is impossible to fit a potential everywhere in continuous space with few point charges; hence more subtle methods should be used. The matrix ||Aatij || is indeed highly singular, in the sense that there are only a few large wk’s, and most are very small; see Figure 2. What is important is that the atoms that are close to the QM system will produce large wk values. Following the SVD theory, we therefore restrict the expansion of the charges to those basis vectors with large singular values L
