Electrostatic and quantum chemical investigation of the proton pumping mechanism of cytochrome c oxidase

ELECTROSTATIC AND QUANTUM CHEMICAL INVESTIGATION OF THE PROTON PUMPING MECHANISM OF CYTOCHROME cOXIDASE DISSERTATION submitted to the Faculty of Biology, Chemistry and Geoscience of the University of Bayreuth, Germany for obtaining the degree of Doctor of Natural Sciences presented by Punnagai Munusami Bayreuth 2008 Die vorliegende Arbeit wurde in dem Zeitraum von Januar 2005 bis Mai 2008 an der Universit¨ at Bayreuth unter der Leitung von Professor G. Matthias Ullmann erstellt. Vollst¨ andiger Abdruck der von der Fakult¨ at Biologie, Chemie und Geowissenschaften der Universit¨ at Bayreuth genehmigten Dissertation zur Erlangung des akademischen Grades Doktor der Naturwissenschaften (Dr. rer. nat.). Erster Pr¨ ufer: Prof. Dr. G. Matthias Ullmann Zweiter Pr¨ ufer: Prof. Dr. Holger Dobbek Dritter Pr¨ ufer: Prof. Dr. Andreas Fery Pr¨ ufungsvorsitz: Prof. Dr. J¨ urgen Senker Tag der Einreichung: 09.06.2008 Kolloquium: 19.09.2008 ELECTROSTATIC AND QUANTUM CHEMICAL INVESTIGATION OF THE PROTON PUMPING MECHANISM OF CYTOCHROME cOXIDASE SUMMARY Cytochrome coxidase is a crucial enzyme in the respiratory chain. It catalyzes the reduction of oxygen to water and utilizes the free energy of the reduction reaction for proton pumping across the inner-mitochondrial membrane, a process which results in a membrane electrochemical proton gradient. For each oxygen molecule, eight protons are taken up from the matrix of the mitochondria. Four protons together with four electrons are required to reduce oxygen to water at the Fea3-CuBbinuclear center and another four protons are translocated across the membrane. Although several high resolution structures have been solved for this enzyme, the molecular mechanism of the proton pumping and electron transfer is not understood. Recent studies on the cytochrome coxidase CuBcenter suggested deprotonation of the CuB bound imidazole ring of histidine (His291 in mammalian cytochrome coxidase or His334 in Rhodobacter sphaeroides cytochrome coxidase) as a key element in the proton pumping mechanism [1–6]. The central feature of this proposed mechanism is that the pKavalue of the imidazole significantly lowered depending on the redox state of the metals in the binuclear center. The energetic feasibility of this mechanism is tested in this work. To find a reliable method to calculate effective charges for the pKacalculations, the charge distribution of the tripeptide Ala-Asn-Ala with different conformations has been analyzed using the hybrid density functional method (B3LYP) with 6-31G* basis set. Population analysis methods (such as Mulliken and Natural Population Analysis) and electrostatic potential (ESP) methods (such as CHELPG, MK and RESP) are used to analyze the charge distribution of the tripeptide Ala-Asn-Ala. This extensive study provided a better understanding of each method and the parameters which influence the partial atomic charges. The results show that ESP methods like CHELPG and MK give reliable charges when proper sampling points are used for the potential fit. To comprehend the role of the CuBbound histidines in the reaction mechanism of cytochrome coxidase, density functional theory is used in combination with continuum electrostatics to calculate the pKavalues of these imidazole rings in the aqueous solution as well as in the protein. The pKavalues of His334, His333 and H2O molecule are calculated both in oxidized and reduced state of CuBcenter. The Finite Difference Poisson Boltzmann (FDPB) method and the conductor-like polarizable continuum model (C-PCM) are used to determine the solvation free energies in aqueous solution. All possible protonation equilibrium reactions in the CuBcenter are studied to understand the deprotonation reactions of the bound H2O molecule, His333 and His334. In aqueous solution, pKavalues of 15.2, 15.9 and 7.4 were obtained for deprotonation of His334, His333 and H2O 7 respectively. These pKavalues in aqueous solution show that His334 and His333 are likely to be protonated at physiological pH. The protein environment shifts the pKavalues of the CuBligands to even higher values in the range between 15 to 60. These pKavalues of CuBligands are significantly higher compared to aqueous solution. The high pKavalues show that His334 is protonated during all steps of the catalytic cycle and demonstrate that the Fe and Cu ion oxidation states do not lower the pKa values of CuBligands and involved in shifting the pKavalues of CuBligands to higher values. These results are incompatible with the proposed role of His334 as a key element in the pumping mechanism. According to the pKavalues, the proton pumping model as suggested by Stuchebrukhov [1] might not be possible with the involvement of His334. The pKavalues of the His333 in the CuBcenter are always shifted to higher values both in the reduced and in the oxidized state of the CuBcenter. The pKavalues of His333 show that this residue is likely to be protonated in the protein and an involvement in the reaction mechanism of cytochrome c oxidase can therefore be ruled out. ZUSAMMENFASSUNG Cytochrom cOxidase ist ein wichtiges Enzym in der Atmungskette. Es katalysiert die Reduktion von Sauerstoff zu Wasser und nutzt die freie Energie der Reduktion, um Protonen durch die innere mitochondriale Membran zu pumpen, ein Vorgang, der zu einem elektrochemischen Protonengradienten ¨ uber der Membran f¨ uhrt. F¨ ur jedes Sauerstoffmolek¨ ul werden acht Protonen von der mitochondrialen Matrix aufgenommen. Vier Protonen zusammen mit vier Elektronen sind n¨ otig, um Sauerstoff am Binuklearzentrum (Heme-Fea3—CuB) zu Wasser zu reduzieren und weitere vier Protonen werden durch die Membran transportiert. Obwohl f¨ ur dieses Enzym einige hochaufgel¨ osten Strukturen bestimmt wurden, ist der molekulare Mechanismus des Protonenpumpens und des Elektronentransfers nicht verstanden. Neuere Studien am CuB-Zentrum der Cytochrom cOxidase legen nahe, dass die Deprotonierung eines am CuBgebundenen Imidazolrings ein Schl¨ usselelement im Protonenpumpmechanismus darstellt (His291 in der S¨ augetier Cytochrom cOxidase oder His334 in der Rhodobacter sphaeroides Cytochrom cOxidase) [1–6]. Der zentrale Punkt dieses vorgeschlagenen Mechanismuses ist eine erhebliche Verschiebung des pKa-Wertes eines Imidazols im CuB-Zentrum zu niedrigeren Werten in Abh¨ angigkeit vom Redoxzustand der Metalle im binuklearen Zentrum. Die energetische M¨ oglichkeit dieses Mechanismus wird in dieser Arbeit gepr¨ uft. Um eine verl¨ assliche Methode f¨ ur die Bestimmung effektiver Ladungen f¨ ur die pKa-Berechnung zu finden, wurde die Ladungsverteilung im Tripeptids Ala-Asn-Ala in verschiedenen Konformationen mittels einer Hybriddichtefunktionsmethode (B3LYP) mit 6-31G* als Basissatz analysiert. Populationsanalysemethoden (Mulliken-Analyse und Natural Population Analysis) und Methoden, die das elektrostatischen Potential (ESP) verwenden (CHELP, MK und RESP), wurden benutzt, um die Ladungsverteilung des Tripeptids Ala-Asn-Ala zu analysieren. Diese ausf¨ uhrlichen Untersuchung der Methoden zur Berechnung der Partialladungen liefert ein besseres Verst¨ andnis jeder Methode und der Parameter, die die partielle Atomladung beeinflussen. Die Ergebnisse zeigen, dass die ESP-Methoden, wie CHELP und MK, verl¨ assliche Ladungen ergeben, wenn geeignete Probenpunkte f¨ ur den Potentialfit benutzt werden. Um die Rolle der an CuBgebundenen Histidine im Reaktionsmechanismus der Cytochrom c Oxidase zu verstehen, wurden Dichtefunktionstheorie-Methoden in Verbindung mit Kontinuumselektrostatik-Rechnungen verwendet, um die pKa-Werte der Imidazolringe in w¨ assriger L¨ osung und auch im Protein zu berechnen. Die pKa-Werte von His334, His333 und des gebundenen Wassermolek¨ uls wurden im oxidierten und reduzierten Zustand des CuBZentrums berechnet. Solvatationsenergien in w¨ assriger L¨ osung wurden mit Hilfe von finite-difference9 16 CONTENTS 3.1.1 Wave function based methods . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.1.2 Potential based methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.2 Structure preparation and Computational Methods . . . . . . . . . . . . . . . . . 72 3.3 ResultsandDiscussion.................................. 75 3.3.1 Mullikencharges.................................. 76 3.3.2 NPAcharges .................................... 76 3.3.3 CHELPGcharges.................................. 76 3.3.4 MKcharges..................................... 76 3.3.5 RESPcharges.................................... 78 3.4 Conclusions......................................... 80 4 Protonation and redox potentials of cytochrome cnitrite reductase 87 4.1 Hemes and the calcium binding site in cytochrome cnitrite reductase . . . . . . 89 4.2 Preparation of the crystal structure of cytochrome cnitrite reductase . . . . . . . 90 4.2.1 Density Functional calculations . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.2.2 Continuum electrostatic calculations . . . . . . . . . . . . . . . . . . . . . . 91 4.3 ResultsandDiscussion.................................. 93 4.4 Conclusions......................................... 99 5 pKacalculations of binuclear center of cytochrome coxidase using DFT calculations 103 5.1 pKacalculations...................................... 104 5.2 ComputationalMethods ................................. 105 5.2.1 Density Functional calculations . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.2.2 Chargefitting.................................... 108 5.2.3 Solvation free energy calculations . . . . . . . . . . . . . . . . . . . . . . . . 109 5.3 ResultsandDiscussion.................................. 111 5.3.1 pKavalues of the CuBcenter........................... 112 5.3.2 pKavalues of the heme a3center ........................ 120 5.4 Conclusions......................................... 121 6 pKacalculations of CuBligands in cytochrome coxidase 125 6.1 Previous computational work on cytochrome coxidase................ 126 6.2 Structure preparation and models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 6.2.1 Preparation of X-ray structure of protein . . . . . . . . . . . . . . . . . . . . 129 6.2.2 Redoxcentermodels................................ 130 CONTENTS 17 6.3 Density Functional calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 6.4 Electrostatic calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.4.1 Calculation of average pKainprotein...................... 135 6.5 ResultsandDiscussion.................................. 137 6.6 Conclusions......................................... 149 7 Concluding Remarks and Outlook 151 Bibliography 153 18 CONTENTS CHAPTER 1 INTRODUCTION 1.1 INTRODUCTION 1.1.1 CYTOCHROME cOXIDASE - A REDOX-DRIVEN MOLECULAR MACHINE Cytochrome coxidase is the terminal member of the electron transport system of mitochondria and many bacteria. It catalyzes the reduction of molecular oxygen to water and pumps protons across the membrane [7, 8]. In this process the membrane electrochemical proton gradient is generated; the energy stored by the proton gradient is subsequently utilized for ATP synthesis [9]. Cytochrome coxidase is responsible for catalyzing the reduction of more than 95% of the oxygen taken up by aerobically growing higher organisms. Cytochrome coxidase is structurally classified as a member of the superfamily of heme-copper containing terminal oxidases composed of 4 to 13 subunits whose largest and most hydrophobic subunits, I, II and III, are encoded by mitochondrial DNA. These subunits collectively have 18 hydrophobic segments forming membrane-spanning helices similar to those of cytochrome b[10]. The protein subunit I of most cytochrome coxidases contain two heme amolecules, called heme aand heme a3and a copper B center (CuB). Heme a3and CuBform a binuclear center where molecular oxygen is reduced into water. Electrons from cytochrome care first transferred to the copper A center (CuA), which is located in the subunit II. The copper A center passes them to heme a, which transfers the electrons to the binuclear center, where the molecular oxygen is reduced into water (see Eq. (1.1)). The reduction of oxygen to water in the catalytic center of cytochrome coxidases generates energy that is necessary for proton translocation from the mitochondrial matrix. 4 cyt cred + 8 H+ (in)+ O2−→ 4 cyt cox + 2 H2O + 4 H+ (out)(1.1) 1.1.2 STRUCTURE AND FUNCTION OF CYTOCHROME cOXIDASE Crystallization and X-ray diffraction analysis of the mitochondrial and of several bacterial cyctochrome coxidases have set the stage for understanding of the complex functions of this enzyme [7]. 19 20 Introduction Q + + + + ATP synthase Complex IV ADP ATP Intermembrane Inner Mitochondrial Membrane H H H O 2H O 2 Complex III 2 Succinate Complex II Succinate Fumarate NADH NAD+ H Complex I NADH Cyt c dehydrogenase dehydrogenase oxidase Cytochrome bc1Cytochrome c space matrix QH2QH2 Fe Fe Cu Cu Cu Cu Q QH Q 2 Q/ Figure 1.1. The mitochondrial electron-transport chain. The electrons are transferred from complex I to complex IV. Electrons are transferred between complexes I and III by the mobile electron carrier coenzyme Q (Q) and from complexes III to VI by the peripheral membrane protein cytochrome c(Cyt c). Complex II also transfers electrons to Q. The pathways of electron transfer (red) and proton pumping (blue) are indicated. CYTOCHROME cOXIDASE IN MITOCHONDRIAL RESPIRATION. The mitochondrion is the site of eukaryotic oxidative metabolism. In oxidative phosphorylation, electrons are transferred from NADH (Nicotinamide Adenine Dinucleotide reduced form) or FADH2(Flavin Adenine Di nucleotide, reduced form) to O2via membrane-bound protein complexes. The free energy of the electron transfer reaction is coupled to ATP synthesis [10]. A schematic representation of electron transport chain is shown in Figure 1.1. Protein complexes embedded in the inner mitochondrial membrane catalyze the electron transfer from NADH to oxygen. These protein complexes are commonly referred as respiratory chain. Some of the compounds are highly mobile (coenzyme Q (Q) and cytochrome c) which shuttle electrons between the transmembrane protein complexes. The free energy released in the redox reactions is stored as electrochemical gradient of protons across the membrane. This proton gradient is utilized to generate ATP. 2 NADH + 2 H++ O2−→ 2 H2O + 2 NAD+(1.2) The respiratory chain in the inner mitochondrial membrane is commonly organized into four transmembrane complexes, namely NADH dehydrogenase (complex I), succinate dehydrogenase (complex II), cytochrome bc1(complex III) and cytochrome coxidase (complex VI). Two mobile electron carriers Q and cytochrome cshuttle electrons between the transmembrane protein complexes. Complex I passes the electrons from NADH to Q. The electron transfer takes place along ironsulfur (FeS) clusters. During the electron transfer, complex I pumps protons from the matrix to the intermembrane space. 1.1. Introduction 21 Complex II also transfers electrons to Q, however in the case of complex II, the source of the electrons is FADH2produced in the citric acid cycle. Since the standard redox potential of FAD is slightly lower than that of Q, complex II does not pump protons from the matrix to the intermembrane space. The electrons from complex I and complex II are shuttled to complex III by reduced Q molecules. Complex III (cytochrome bc1) passes the electrons to the next mobile carrier cytochrome c which shuttles the electrons to complex IV (cytochrome coxidase). Complex IV passes the electrons to diatomic oxygen and thereby reducing it into two water molecules (final electron acceptor). Both cytochrome bc1and cytochrome coxidase translocate protons from the matrix to the intermembrane space. THE THREE-DIMENSIONAL STRUCTURE OF HEME-COPPER OXIDASES Heme-copper oxidases are membrane proteins found in the respiratory chain of aerobic organisms. They are the terminal electron acceptors coupling the translocation of protons across the membrane with the reduction of oxygen to water. The heme-copper oxidase superfamily is divided into two main branches based on the identities of the electron donating substrates: cytochrome coxidases use a water soluble protein, cytochrome c, as electron donor whereas ubiquinol oxidases use a membrane soluble ubiquinol molecule as their electron donor. Both the members of the superfamily share important structural and functional features. STRUCTURE OF THE CYTOCHROME cOXIDASE FROM Paracoccus denitrificans. The protein was originally isolated as a two-subunit enzyme complex from the cytoplasmic membrane of the soil bacterium P. denitrificans [11]. Michel and co-workers determined the structures of the reduced and oxidized P. denitrificans enzyme and found no structural difference [12]. The four subunits of oxidase has been crystallized in the presence of dodecyl maltoside as a complex with a monoclonal antibody fragment (Fv) directed against an epitope on the hydrophilic domain of subunit II and its structure was determined at 2.8 ˚ A [12, 13]. Later, the two-subunit complex structure was solved at 2.7 ˚ A again using the Fvapproach to increase the polar surfaces of the protein complex and undecyl maltoside as detergent [14]. STRUCTURE OF THE CYTOCHROME cOXIDASE FROM Rhodobacter sphaeroides. The structure of cytochrome coxidase from R. sphaeroides has been solved by Iwata et al. [15]. The crystal structures were determined for the wild type and a mutant by replacing the glutamate residue 286 of subunit I by glutamine (see Figure 1.2). STRUCTURE OF THE CYTOCHROME cOXIDASE FROM BOVINE HEART. The crystal structure of bovine heart cytochrome coxidase at 2.8 ˚ A was determined by Yoshikawa et al. [16–18] in 1995. The crystal structures of the bovine heart cytochrome coxidase in the fully oxidized and the fully reduced states were determined by the same group in 1998 [16]. The protein is composed of 13 different subunits. UBIQUINOL OXIDASE. Ubiquinol oxidases take part in the bacterial electron transport chain oxidizing ubiquinol into ubiquinone and reducing oxygen to water. These enzymes are one 22 Introduction set of the many alternative terminal oxidases in the branched prokaryotic electron transport chain. The overall structure of the ubiquinol oxidase is similar to that of the mammalian cytochrome coxidase, with the addition of a polar ubiquinol-binding site embedded in the membrane. The cytochrome coxidase contains four redox centers: CuA, heme a, heme a3and CuB. The oxygen reduction takes place in the binuclear center Fea3-CuBand utilizes cytochrome cas an electron donor. A second major branch of the aa3-cytochrome coxidase uses ubiquinol or menaquinol as the reducing substrate, and in these cases the CuAcenter is absent. Such enzymes are the cytochrome bo3of Escherichia coli [19, 20] and cytochrome aa3-600 of Bacillus subtilis [21–23]. The cytochrome aa3-type terminal quinol oxidase of B. subtilis catalyzes the four-electron reduction of oxygen into water. STRUCTURE OF THE ba3-CYTOCHROME cOXIDASE FROM Thermus thermophilus. The structure of the ba3-cytochrome coxidase from T. thermophilus has been reported at a resolution of 2.4 ˚ A [24]. The crystal structures of recombinant cytochrome ba3-cytochrome coxidase from T. thermophilus was reported at 2.3 ˚ A by Hunsicker-Wang and co-workers [25]. The model of the ba3-cytochrome coxidase is composed of three protein subunits I, II and IIa. The main part of the complex is formed by subunit I with 13 transmembrane helices, which binds the heme b, heme as3as well as CuB. The a-type heme of the ba3-oxidase corresponds to the heme as present in the SoxB-type terminal oxidases [24]. Heme bis the simplest protoheme containing a low-spin iron with two histidine residues as axial ligands. The binuclear center is formed between heme as3and CuBas in the cytochrome coxidase. STRUCTURE OF THE bo3-CYTOCHROME cOXIDASE FROM Escherichia coli. The structure of cytochrome bo3oxidase from E. coli was reported at 3.5 ˚ A by Abramson et al. [26]. The cytochrome bo3ubiquinol oxidase is a four-subunit heme-copper oxidase that catalyzes the four-electron reduction of O2to water and functions as a proton pump. All the redox centers are located in subunit I, with a low spin protoheme bacting as an electron donor to a binuclear center that is composed of an o-type heme, heme o3, and a copper ion CuB. Subunits of I, II and III of ubiquinol oxidase are homologous to the corresponding subunits in the aa3-type cytochrome coxidase, and the ligands of the two heme groups and the CuBhave been identified as histidine residues. In contrast to the cytochrome coxidase, subunit II of the ubiquinol oxidase has neither CuAcenter, nor a cytochrome cbinding site. Instead, the heme breceives electrons directly from a membrane solubilized ubiquinol molecule. SUBUNITS OF CYTOCHROME cOXIDASE The subunits of cytochrome coxidase from R. sphaeroides are discussed in the following section. Subunit I (see Figure 1.2) of bacterial cytochrome coxidase is largely embedded in the membrane, with its 12 transmembrane helices shaped in a three-winged propeller arrangement [7, 12]. The N-terminus and the long, exposed C-terminus of the polypeptide face the cytoplasmic side. The three redox centers, the two a-type hemes and the copper B center, are ligated by amino acid side chains of this subunit. Histidines are the axial ligands to the lowspin heme a, whereas a histidine and a presumed hydroxyl or a water molecule are ligands 1.1. Introduction 23 2 H O O2 SU II SU III SU I + CuA heme a 3CuB 4H + 4H 4H + Cyc c4 matrix membrane Intermembrane space SU IV heme a Figure 1.2. The X-ray structure of cytochrome coxidase from Rhodobacter sphaeroides.The electrons are transferred from cytochrome cto copper A center. Protons are pumped from matrix to intermembrane space and four more protons are delivered to the binuclear center. Oxygen is reduced to water in the binuclear center (heme a3...CuB). to the high-spin heme a3moiety. Both hemes are oriented perpendicular to the membrane plane. Heme a3, together with a copper ion (CuB)in its immediate vicinity, forms the binuclear center where oxygen binding and reduction takes place. Subunit II (see Figure 1.2) of bacterial cytochrome coxidase has a bipartite structure. The N-terminal has two transmembrane helices followed by a hydrophilic, 10-stranded β-barrel domain extending into the periplasm, comprising the CuAcenter. The CuAcenter contains two copper ions in a mixed-valence (CuI.CuII)state and 2.6 ˚ A apart, giving rise to a characteristic EPR (electron paramagnetic resonance) signal in the oxidized state which was observed in the other enzymes as well [27]. The two copper ions are bridged by two cystein thiolates. Subunit III (see Figure 1.2) of bacterial cytochrome coxidase is fully embedded in the membrane. No redox cofactors are associated with this subunit. When both subunits III and IV are removed from the P. denitrificans cytochrome coxidase there is no loss in the catalytic functions. It has been suggested that subunit III stabilizes the integrity of the binuclear center in subunit I [28]. The cleft of the subunit III may be a binding site for other membrane proteins. Subunit III may be involved in assembly of the oxidase or form the entrance to a oxygen channel leading to the active site [29]. It possesses seven transmembrane helices that are divided by a large V-shaped cleft into two bundles, one formed by the first two helices, and the other by helices III to VII. In this cleft, lipid molecules are found to be firmly bound to the 24 Introduction conserved residues. In mammalian oxidase three hydrophobic channels were proposed and these hydrophobic channels are suggested to be the potential pathways for oxygen to reach the binuclear center [18]. The oxygen channels start at the protein-membrane interface near the center of the lipid bilayer, where oxygen solubility is much higher than in the aqueous phase. One of such channel has also been identified in the bacterial oxidase [29]. It starts in the V-shaped cleft of the subunit III directly above a tightly bound lipid molecule and leads through subunit I into the binuclear site. Subunit IV (see Figure 1.2) of the bacterial enzyme of both P. denitrificans and R. sphaeroides consists of a single transmembrane helix in contact with subunit I and III. The function of this small subunit is unknown [12]. THE REDOX CENTERS OF CYTOCHROME cOXIDASE C256 H217 H260 E254 CuA M263 C252 Figure 1.3. The copper A (CuA) center of cytochrome coxidase from R. sphaeroides. The copper A (CuA) center with its ligating residues is shown. The electrons are transferred from cytochrome cto the copper A center. The amino acid numbering of R. sphaeroides is used throughout this thesis. COPPER A(CuA)CENTER. The CuAcenter is located 8 ˚ A above the membrane surface. The CuAcenter contains two copper ions (see Figure 1.3). These copper ions are bridged by two cysteins sulfur atoms and have additional protein ligands. The two copper ions of the CuA center are coordinated by two His, one Met, a backbone carbonyl oxygen of a Glu and two bridging Cys residues. The spectroscopic measurements indicate that in the reduced form of the CuAcenter both copper ions are in their Cu(I) state whereas in the fully oxidized form, the newly acquired electron appears to be delocalized between the two copper ions such that they assume the mixed [Cu1.5+...Cu1.5+]state. 1.1. Introduction 25 H102 heme a H421 Figure 1.4. The heme aof cytochrome coxidase from R. sphaeroides.The heme a porphyrin ring system is shown with two histidine residues as axial Fe-ligands. The amino acid numbering of R. sphaeroides is used throughout this thesis. HEME aCENTER. The heme acenter consists of two histidine residues as axial iron-ligands (see Figure 1.4). It transfers the electrons from CuAto the binuclear center. The heme is non-covalently bound to the protein and heme acontains a formyl group and a hydrophobic hydroxyethyl-farnesyl group. The heme environments in P. denitrificans and bovine heart cytochrome coxidase are very similar. HEME a3AND COPPER B(CuB)CENTER: THE BINUCLEAR CENTER. The binuclear center is formed by heme a3and CuBwhich is the catalytic center for O2reduction (see Figure 1.5). The heme a3iron appears to be five fold coordinated with a histidine of subunit I. The copper ion in the CuBcenter is coordinated by Natoms of His334 and His333 and Nδatom of His284 and the histidines are arranged in an equilateral triangle, centered on CuB. Molecular oxygen is supposed to bind between the heme a3iron and CuB. The iron of heme a3is 0.36 ˚ A out of the heme plane in P. denitrificans but almost within the plane in the bovine heart. Crystallographic studies on cytochrome coxidase have revealed an unique and unexpected posttranslational modification in the enzyme active site [16]. A tyrosine (Y288) which is located very near to the binuclear center is covalently linked to the Nof histidine residue (H284) which also serves as a ligand to the CuBcenter (see Figure 1.5). The cross-linked tyrosine and histidine was observed in both bacterial and mammalian cytochrome coxidases. AMINO ACIDS AS REDOX CENTERS. In additional to four metal redox centers, there are other amino acids near the active site of the enzyme that could conceivably form radicals and act as a redox-active centers in cytochrome coxidase. The tyrosine (Y288) is the most proximate amino acid which is cross-linked to the CuBhistidine (H284) [30]. The EPR evidence 32 Introduction ASTATE. The proposed catalytic cycle for the molecular steps taking place in the active site of cytochrome coxidase during catalysis is given in Figure 1.8. Oxygen binds to the reduced binuclear center and forming the A state (adduct after oxygen binding). PSTATE. A state leads to the so-called peroxy state P. The P state appears to exist in two forms, Pm, a two electron reduced and the Prstate. The Pmstate resulted from the internal electron redistribution and the tyrosine is proposed to be neutral radical (see Figure 1.8). Weng and Baker [55] were the first to suggest that the P state is not a peroxy state but the oxo-ferryl state with hydroxyl group bound to the CuBcenter. They suggested that a tryptophan is the source of the missing electron in analogy to cytochrome cperoxidase. Later Kitagawa and co-workers [56] provided evidence by Raman spectroscopy that P state is a hydrogen-bonded oxo-ferryl state. The structure of the Pmstate shown in Figure 1.8 has been proposed several times [57]. The existence of the covalent tyrosine-histidine cross-link should be taken as the evidence that a tyrosine radical is formed during the catalytic cycle of cytochrome coxidase, because the cross-linking of tyrosines is typical of a radical reaction catalyzed by peroxidase. The Prstate of the binuclear center may be considered as a “highenergy” state, the decomposition of which will drive proton translocation [58]. The subsequent transformation of the Prstate into F state is coupled to translocation of the proton across the membrane, net uptake of another proton into the binuclear center and an accompanying generation of transmembrane electric potential [59]. FSTATE TO OSTATE. The F state results from the third electron transfer from cytochrome c together with the acquisition of two protons which convert the tyrosine radical to phenolate state. There are overwhelming spectroscopic data that heme a3is in the oxyferryl form in F state (see Figure 1.8). However, there appears to be multiple forms and ambiguity with respect to the redox status and the protonation state of groups in the immediate vicinity of the active site, depending on how F state is generated. A fourth and final electron transfer and proton acquisition yields the oxidized O state via H state. The H state intermediate possesses hydroxyl group at the heme a3. OAND RSTATE. In O state, all the four redox centers are fully oxidized and the heme a3iron is bound with the water ligand and CuBcenter is bound with hydroxyl group. The tyrosine is covalently cross-linked to the CuBligand of His284. The oxidized binuclear complex is reduced to R state via the formation of the one-electron reduced E state. Protons are translocated from matrix during this process. Wikstr¨ om and co-workers worked to establish much of the conceptual and an experimental framework for studying the oxidase mechanism, including important experimental evidence that proton pump is coupled only to the P→F and F→O state transitions and that each of these one electron redox steps results in the pumping of two protons [8, 60]. This paradigm has been generally accepted for the past decade, though some problems have been pointed out and discussed [61]. Now the paradigm that pumping is coupled only to the last two steps in the catalytic cycle is being challenged by Michel [53] on the basis of a re-evaluation of several key experiments. The conclusion of the re-evaluation is that the F →O state transition is coupled to the pumping of only one proton and not two. Although considerable progress has 1.3. Outline of the thesis 33 been made concerning the proton input pathways in cytochrome coxidase, very little is known about the mechanism and how the proton pump actually works. PROTON–COUPLED ELECTRON TRANSFER The cytochrome oxidase energetically couple the electron transfer reactions associated with reduction of oxygen to water and pumps proton across the membrane. Even though a vast amount of structural and functional information of cytochrome coxidase are available from experimental and theoretical data [8, 14, 16, 50–52], actual step of coupling the the redox reactions to the proton translocation is poorly understood. How the protons are delivered exactly to the binuclear center, the proton translocation and exit pathways and the chemical intermediates involved in the catalytic cycle are still remain unclear. 1.3 OUTLINE OF THE THESIS The aim of the thesis is to understand the proton pumping mechanism of cytochrome coxidase. Although the structures of cytochrome coxidase has been solved for several organism, the molecular mechanism of proton pumping remains unclear. In this thesis, the reaction mechanism of cytochrome coxidase is analyzed by combining the density functional theory (DFT) and continuum electrostatic calculations. The theory behind the electrostatic calculations, Poisson-Boltzmann equation, titration behavior calculations and DFT are described in chapter 2 and different charge methods are reported in chapter 3. Protonation probabilities and redox potentials of cytochrome cnitrite reductase which serves as a simpler electron transfer system to study cytochrome coxidase are discussed in chapter 4. To investigate the reaction mechanism and the role of the histidines bound with the CuB center and the histidine coordinating the heme a3, the pKacalculations were performed on the CuBand heme a3using DFT in combination with continuum electrostatic models. The pKavalues are calculated using different basis sets and different solvation models. The Finite Difference Poisson Boltzmann (FDPB) method and the conductor-like polarizable continuum model (C-PCM) are used to determine the solvation free energies in aqueous solution. The influence of different charges, basis sets and solvation models on the pKacalculations of the CuBhave been studied (Chapter 5). The pKavalues of the CuBligands and heme a3center in aqueous solution are reported in Chapter 5. In chapter 6 the average pKavalues of the CuBligands in the presence of different heme a3 redox states in cytochrome coxidase are reported. This thesis provides insights into the role of CuBligands in the reaction mechanism of cytochrome coxidase. The most accurate methods are reported for the pKacalculations of the CuBand heme a3center. This thesis is thus a step towards a better understanding of the role of the CuBligand in the reaction mechanism of cytochrome coxidase. 34 Introduction CHAPTER 2 COMPUTATIONAL METHODS 2.1 ACID-BASE AND REDOX REACTIONS EQUILIBRIA 2.1.1 FUNDAMENTAL DESCRIPTION OF ACID-BASE AND REDOX EQUILIBRIA Biological molecules, such as proteins and nucleic acids bear numerous functional groups such as carboxyl and amino groups that can undergo acid-base reactions. The dissociation of a proton from a monoprotic acid is generally given by: HA A−+ H+(2.1) The free energy change (∆Ga) of this reaction can be related to the equilibrium constant (Ka), ∆Ga=−RT ln Ka(2.2) Ka=[H+][A−] [HA] (2.3) The pH of the solution is defined as the negative decadic logarithm of the hydrogen ion concentration and the pKaof an acid is defined as the negative decadic logarithm of the Ka values. pH = −log[H+](2.4) pKa=−logKa(2.5) The Henderson-Hasselbalch equation combines Eq. (2.4) and (2.5) to pH = pKa+ log [A−] [HA] (2.6) The protonation probability is given by [62]: pprot =[HA] [HA] + [A−]=10pKa−pH 1 + 10pKa−pH (2.7) 35 36 Computational methods The pKavalue of an acid is the pH value at which the concentration of the protonated and deprotonated forms of the acid equal. Similar to the above equilibrium protonation reaction, the equilibrium between the redox reaction is, Aox +e−A− red (2.8) The reduction equilibrium constant KET for this reaction is, KET =[A− red] [Aox][e−](2.9) Eq. (2.9) is analogous Eq. (2.6). The standard redox potential of the above redox reaction is given by: E0=RT Fln KET (2.10) The redox potential of the solution is given by, E=−RT Fln[e−](2.11) The Nernst equation combines the standard redox potential E0and the redox potential of the solution E E=E0+RT Fln [A] [A−](2.12) where Fis the Faraday constant, Ris the gas constant and Tis the temperature. The probability of finding group A in the reduced state is therefore given by [62]: pred =[A−] [A−] + [A] =exp RT F(E0−E) 1 + exp RT F(E0−E)(2.13) 2.1.2 COMPUTATION OF ACID-BASE AND REDOX EQUILIBRIA The pKais directly related to the free energy of the deprotonation reaction in aqueous solution ∆Gdepro water by the following Eq. (2.14): pKa=1 ln 10kBT∆Gdepro water (2.14) The ∆Gdepro water can be expressed as a sum of two contributions: the solvation energy difference ∆∆Gdepro solv between the associated and the dissociated system and the gas phase deprotonation energy ∆Gdepro vac . These terms can be obtained from the thermodynamic cycle (see Figure 2.1). pKa=1 ln 10kBT(∆Gdepro vac + ∆∆Gdepro solv )(2.15) 2.1. Acid-base and redox reactions equilibria 37 (A ) (H ) + water vacuum AH AH AH+ ∆G ∆G∆G∆G(AH) solv depro vac solv solv AH+ G∆depro water Figure 2.1. Thermodynamic cycle to calculate absolute pKavalues. The free energy of dissociation is calculated in vacuum and the reactant and products are then transferred from vacuum to water. The free energy of dissociation of a proton from an acid in water (∆Gdepro solv ) is calculated indirectly. The solvation energy difference ∆∆Gdepro solv is obtained from Eq. (2.16) ∆∆Gdepro solv = ∆Gsolv(A−)+∆Gsolv(H+)−∆Gsolv(AH)(2.16) The solvation energy of the protonated and deprotonated states, ∆Gsolv(A−)and ∆Gsolv(AH) can be calculated by solving the Poisson-Boltzmann equation Eq. (2.19). The solvation energy of a proton is measured experimentally from the potential of the standard hydrogen electrode. The solvation energy of –264.6 kcal/mol [63] is used for the proton in the present study. The gas phase protonation energy ∆Gdepro vac is given by Eq. (2.17) ∆Gdepro vac = ∆Hdepro vac + ∆Hdepro vib +Htrans(H+) + ∆(pV )−T[S(H+)] (2.17) where ∆Hdeprot vac is the difference in the gas phase energy of the associated (protonated) and dissociated (deprotonated and hydrogen ion) system which can be obtained from quantum chemical calculations. The ∆Hdeprot vib is the change in the vibrational energy between the protonated and deprotonated states and this can be obtained from normal mode analysis. The Htrans(H+) = 3 2RT is the translation energy of a proton and ∆(pV )is the energy change due to the volume change in the gas phase reaction which is estimated to be kBTfrom the ideal gas approximation. The T[S(H+)] is the entropic contribution of the gas phase free energy of proton and is set to 7.8 kcal/mol [62] as derived from the Sackur-Tetrode equation. The redox potential Eredox 0can also be computed by a similar approach. The redox potential Eredox 0can be calculated from Eq. (2.18) Eredox 0=1 F(∆Hredox vac + ∆∆Gredox solv )+∆SHE (2.18) 38 Computational methods where ∆Hredox vac is the difference in the gas phase energy between the oxidized and reduced states. The ∆∆redox solv is the solvation energy difference between the oxidized and reduced states which can be calculated by solving the Poisson-Boltzmann equation. Fis the Faraday constant (23.06 kcal/mol) and ∆SHE is the standard potential of the hydrogen electrode (–4.43 V) [62]. 2.2 CONTINUUM ELECTROSTATICS Electrostatic effects are very important in many phenomena in biology and chemistry. In biology, electrostatics is a key component in determining the structure of proteins, nucleic acids and membranes. Electrostatics of proteins can be studied by different models. The quantum chemical calculations can give accurate results, but it is computationally expensive and time consuming method. An alternative approach involves continuum or macroscopic models, in which the solvent properties are described in terms of average values. By this model it is possible to describe solute molecules in atomic details and treating the solvent in terms of average properties [64]. Following approximations will be used in the continuum electrostatics model. The solvent water is assumed as a continuum liquid with high dielectric constant (= 80) (see Figure 2.2). The protein system will be assigned a low dielectric constant (= 2 −4). The solvent which contains mobile ions are modeled as mobile charges. The protein system is represented as a low dielectric medium containing fixed point charges and the mobile ions are assumed to be infinitely small in size [62, 65, 66] (see Figure 2.2). + + − + + + + − − − − + ++ + + − − − − − mobile ions Solvent with with ionic strength I Fixed charges within the molecule (low ) ε (high )ε Ion exclusion layer Figure 2.2. A molecule in a heterogeneous dielectric medium. A molecule of low dielectric with fixed charges in a high dielectric solvent with mobile charges (ion). 2.2. Continuum electrostatics 39 2.2.1 THE POISSON-BOLTZMANN EQUATION The treatment of charges in vacuum is the simplest form of the classical electrostatics. Poisson equation can be applied to treat such electrostatic interactions. The electrostatic potential φ(~r)at position (~r), arising from charge density ρ(~r)in vacuum is given by, ~ ∇2φ(~r) = −ρ(~r) ε0 (2.19) where ε0is the dielectric constant in vacuum and ~r is given in Cartesian coordinates. The ∇is the Nabla-operator (~ ∇= (∂/∂x, ∂/∂y, ∂/∂z)) and ~ ∇2=~ ∇ · ~ ∇is Laplace-operator. The Solutions of Eq. (2.19) yields the familiar coulombs law. For inhomogeneous dielectric medium, the dielectrium varies with position (~r), and Poisson equation is modified to ~ ∇ · [ε(~r)~ ∇φ(~r)] = −ρ(~r) ε0 (2.20) The protein system is assumed to have mobile ions. The charge distribution of the mobile ions ρions(~r)can be described by Debye-H¨ uckel theory. The theory describes the distribution of ions by combining electrostatics with statistical thermodynamics and uses the Boltzmann statistics based on the electrostatic interaction energy Wibetween ion iand the solute [67]. The charge distribution resulting from the distribution of ions ci(~r)and other ions are ρions(~r) = K X i=1 ci(~r)·Zi·e0= K X i=1 cbulk iZie0·exp −Wi(~r) RT (2.21) where cbulk iis the concentration of ion iin the bulk, Kis the number of different type of ions, Ziis the unitless formal charge of ion type i,e0is the elementary charge, Ris the universal gas constant and Tis the absolute temperature. The summation runs over all Ktypes of ion i. The interaction energy Wiis approximated by the potential energy Epot of ion iin the electrostatic potential φ(~r)due to charge distributions of solute ρsolute(~r)and mobile ions ρions(~r). This approximation Wi≈Epot =Zie·φ(~r)yields ρions(~r) = K X i=1 cbulk iZie0·exp −Zie0·φ(~r) RT (2.22) With ρ(~r) = ρsolute(~r) + ρions(~r), Eq. (2.19) can be re-written to give the Poisson-Boltzmann equation (PBE) ~ ∇ · [ε(~r)~ ∇φ(~r)] = 1 ε0 ρsolute(~r) + K X i cbulk iZie0·exp −Zie0·φ(~r) RT !(2.23) 40 Computational methods Using the linear approximation exp(−x)≈1−xvalid for small x1, for Zie0·φ(~r)/RT 1 K X i cbulk iZie0·exp −Zie0·φ(~r) RT = K X i cbulk iZie0− K X i cbulk iZie0·Zie0·φ(~r) RT (2.24) As the total charge of the ions in solution is equal to zero, the first term on the right side of Eq. (2.24) is zero. ~ ∇ · [ε(~r)~ ∇φ(~r)] = 1 ε0 ρsolute(~r) + K X icbulk iZ2 ie2 0 RT φ(~r)!(2.25) The linear Poisson-Boltzmann equation (LPBE) can be rewritten by introducing the ionic strength and the Debye-H¨ uckel inverse, κ I=1 2 K X i=1 cbulk iZ2 i(2.26) κ2(~r) = 8πe2 0I(~r) RT (2.27) ~ ∇ · [ε(~r)~ ∇φ(~r)] = −4πρsolute(~r) + κ2φ(~r)(2.28) Eq. (2.28) is solved numerically to yield φ(~r). The Poisson-Boltzmann equation is quite good for the mobile ions (monovalent and multivalent). The Linear Poisson-Boltzmann equation breakdowns for solvent with high concentration but these limitations will not affect results for the biological system since both ionic strengths and the electrostatic interactions are moderate in these systems. NUMERICAL SOLUTION TO LINEAR POISSON-BOLTZMANN EQUATION. Analytical solutions are known only for simple geometries. For complex systems, numerical methods are applied to solve the Poisson-Boltzmann equation. A wide variety of problems have been studied using the finite difference method. However there are also other numerical methods such as boundary element method which tessellates the dielectric boundary between the interior and exterior of the molecule. A finite difference method is applied in the present work to solve the PBE. In finite difference method (see Figure 2.3), the protein (solute) is mapped into a cubic grid along with the solvent. The values of the electrostatic potential, charge density, grid constant h, dielectric constant and ionic strength are assigned to each grid point [68]. The atomic charges usually do not coincide with the grid points and so the charge is allocated to the eight surrounding grid points in such a way that the closer the charges to the grid point, the greater the proportion of its total charge that is allocated. The linear Poisson-Boltzmann 2.2. Continuum electrostatics 41 φ0 φ3 5 ε 1 ε6 φ2φ4 0 0 1 1 0 0 0 1 1 1 0 0 1 1 0 0 0 1 1 1 0 0 1 1 00 00 11 11 0 0 1 1 h φ5 φ1 φ6 ε2ε ε3 ε4 q0 0 Ι Figure 2.3. The cube used in the finite difference method for solving the PoissonBoltzmann equation. equation Eq. (2.28) is integrated over the volume Vof the cubic grid: Z~ ∇ · [ε(~r)~ ∇φ(~r)]dr −Zκ2 0(~r)φ(~r)+4πZρsolute(~r)d~r = 0 (2.29) In the first step, the integral is transformed into a surface integral using Gauss’s theorem: The volume integral is given by: Z~ ∇ · [ε(~r)~ ∇φ(~r)]dA −h3κ2 0φ0+ 4πq0= 0 (2.30) The surface integral is calculated separately for all six sides of the cubic grid element. Z~ ∇ · [ε(~r)~ ∇φ(~r)]dA = 6 X i=1 h2εi(φi−φ0) h(2.31) 6 X i=1 hεi(φi−φ0)−h3κ2 0φ0+ 4πq0= 0 (2.32) the finite difference expression for φ0: φ0=(P6 i=1 hεiφi)+4πq0 P6 i=1 hεi+h3κ2 0 (2.33) 48 Computational methods pKp r=1 RTln10(Go p−Go r)(2.45) The macroscopic p ¯ Kkvalues characterize the equilibria between two macroscopic protonation states of the protein i.e., between the two macrostates with kand k−1protons bound. The kth macroscopic p ¯ Kkvalue is given by p¯ Kk=−log PM mP2N xδ(k−1)e−βGo x,m PM mP2N xδ(k)e−βGo x,m !(2.46) A macroscopic p ¯ Kkvalue is the thermodynamic average over all microscopic equilibria involved in the release of the kth proton. The product over the first to the kth macroscopic ¯ Kk values equals to the kth polynomial coefficient PN mP2N xδ(k)e−βGo x,m of the partition function in Eq. (2.42). The macroscopic p ¯ Kkvalues can consequently also be used to formulate the partition function. Moreover, the total titration curve of a polyprotic acid is often interpreted in terms of macroscopic p ¯ Kkvalues. The pKavalues for three proton binding sites are also calculated in similar way. The system with three proton binding sites are shown in Figure 2.7 [75]. CALCULATION OF TITRATION CURVES The protonation state of a protein with Ntitratable group is 2N. The titration curve of a single site in a protein is given by the thermodynamic average over all possible protonation states at each pH-value hxii= 2N X n=0 x(n) iexp −G(n)(pH) RT  2N X n=0 exp −G(n)(pH) RT  (2.47) To obtain the protonation state energy Gn, the Eq. (2.47) should be evaluated for 2Ntimes. Since cytochrome coxidase and other model proteins are huge, the thermodynamic average hxiigrows exponentially with the number of sites. The average behavior of such system can be obtained by Metropolis Monte Carlo approach which is used to determine the protonation of many interacting sites as a function of pH. THE METROPOLIS MONTE CARLO APPROACH. The Metropolis Monte Carlo method is used to obtain the titration curves of amino acids within the protein [76]. Initially a random protonation state nis chosen, the energy of the corresponding protonation state is calculated. Each Monte Carlo (MC) step chooses the protonation state randomly and the energy of the current protonation state is calculated and the change in energy is computed, i.e., the energy difference between current protonation state and the previous protonation state ∆G = Gold−Gcurrent. 2.2. Continuum electrostatics 49 The current state is accepted with a probability paccording to Metropolis criterion p=   1 if ∆G≤0, exp (−∆G RT ) if ∆G > 0(2.48) The above process is repeated for each MC step. The collection of protonation states start after the equilibration step. The steps involved in MC are shown in Figure 2.8. The Metropolis Monte Carlo approach also effectively sample the protonation states of the strongly coupled sites. The strongly coupled groups have the following states: the singly protonated state (1,0) and (0,1) the low energy states and state (0,0) and (1,1) high energy states. Therefore two transitions are needed to change from (1,0) to (0,1) to avoid the highenergy intermediate states (0,0) and (1,1). To solve this sampling problem, double and triple moves are introduced, which simultaneously change the protonation states of two or three sites in MC step. CORRELATION CURVES. The coupling between the protonation form of titratable group iand group jis given by: c=hxi,ji−hxii·hxji(2.49) where the hxiiand hxjiare the probability of the titratable group iand group jwhich are protonated respectively and hxi,j iis the probability for both are protonated in the same time. The correlation between two sites can be obtained from above Eq. (2.49). TANFORD-ROXBY PKaVALUES The Tanford-Roxby approximation assumes that the average protonation of a titratable residue depends on the average charges of all other titratable groups. The average pKavalue of the residue iinside the protein [77, 78] hpKintr a,i iprot is obtained from Eq (2.50). pKa,i =pKintr a,i + N X j=1 (hxji − x(0) j)Wij (2.50) where the intrinsic pKavalue (pKintr a,i ) is the pKavalue that the particular titratable group would have, if all other titratable groups are in their reference form. This term includes the solvation energy and the interaction with non-titrating residues and the protein backbone. The term Wij represents the interaction energy between the titratable groups iand jin their charged from; hxjirepresents the protonation probability of the group jwhich is obtained by a thermodynamic average over all possible protonation states and x0 jis the reference protonation form of site j. The hxjiare obtained from MC calculations. These average pKavalues do not represent an equilibrium situation, but they are a good approximation of the real pKa value of the site. 50 Computational methods ∆ G= G − G ∆G State Choose random Initial Protonation Protonation of ∆old Change Check 0 Choose random r1 0r Accept new State ∆G Yes No ∆ Keep old State Save State More MC Calculate average Probability Exit Yes No No Yes Site i Site i new Check er /k T B Figure 2.8. The Metropolis Monte Carlo (MC) approach. Initial random protonation state iis chosen and the protonation is changed. The energy difference between the new state and old state is calculated. The new state is accepted if the energy is smaller than zero. If the energy is greater than zero, the Boltzmann factor of the change in energy is checked. If the change in energy is larger than a number randomly chosen between 0 and 1 then the state is accepted. The steps are repeated until convergence of the sampled properties is achieved. 2.3. Quantum Chemistry 51 2.3 QUANTUM CHEMISTRY The application of quantum-mechanical principles to chemical problems has gained momentum in past few decades. The quantum mechanical methods can be of great use in the simulation of proteins. They can be used to obtain parameters for molecular mechanics type potential energy functions. In the present work, the quantum chemical (QC) methods are used to calculate the free energies of the reactions, vibrational energies, solvation energies and partial atomic charges. This chapter described some of the basic principles, methods of quantum chemistry that lead to understand the electronic structure behavior. SCHR ¨ ODINGER EQUATION The principles of density functional theory are conveniently expounded by making reference to conventional wave function theory. The time independent, non-relativistic Schr¨ odinger equation for a many electron system can be written as: HΨ = EΨ(2.51) His the Hamiltonian operator for a system of nuclei and electrons described by position vectors ~ RAand ~ri, respectively. Ψis the wave function of the many particle system and E is the energy of the eigenfunction. In SI units the Schr¨ odinger equation of the hydrogen atom can be written as: −~2 2m∇2 iΨi(~r)−e2 4πε0rΨi(~r) = EiΨi(~r)(2.52) In Eq. (2.52), first term in the left hand side (LHS) is the kinetic energy of the electron system and the second term is the potential energy term and Eiis the energy of the system, where eand mare the electron charge and mass respectively, ~is Plank’s constant divided by 2π, ris the distance between electron and proton and ε0is the permittivity of free space. In the SI system of units the above equation is rather unwieldy to deal with because of the tiny numbers associated with quantities as m,e2and ~2. The inconvenience can be removed by rewriting the equation in terms of ‘natural’ atomic units. The Schr¨ odinger equation can be expressed in atomic units where me=|e|=~= 4πε0= 1 and Eq. (2.52) becomes Eq. (2.53) −1 2∇2 i−1 rΨi(~r) = EiΨi(~r)(2.53) 52 Computational methods The Hamiltonian for Nelectrons and Mnuclei can be expressed as H=− N X i=1 −1 2∇2 i− M X A=1 1 2MA ∇2 A− N X i=1 M X A=1 ZA riA + N X i=1 N X j>i 1 rij + N X A=1 N X B>A ZAZB RAB (2.54) In Eq. (2.54) the distance between the i-th electron and the A-th nucleus is ~riA=|~ri-~ RiA|, the distance between the i-th and j-th electron is ~rij=|~ri-~rj|and the distance between the A-th and B-th nucleus is ~ Rij=|~ Ri-~ Rj|.MAis the ratio of the mass of nucleus Ato the mass of an electron and ZAis the atomic number of nucleus A. The Laplacian operators ∇2 iand ∇2 A involve differentiation with respect to the coordinates of the i-th electron and A-th nucleus, respectively. The first term in the above equation is the operator for the kinetic energy of the electrons; the second is the operator for the kinetic energy of the nucleus; the third term represents the Coulomb attraction between electrons and nuclei; the fourth and fifth term represents the repulsion between electrons and between nuclei, respectively. The Born-Oppenheimer approximation [79, 80] considers electrons as moving particles in the field of fixed nuclei. Due to significant differences in mass between electrons (mel∼10−31 kg) and the nucleus (mnuc∼10−27 kg), the electronic degrees of freedom can be considered to respond instantaneously to any change in the nuclear configuration i.e., their wave functions always correspond to a stationary state. Hence, the interaction between different stationary states. Application of the Born-Oppenheimer approximation enables one to formulate the electronic Hamiltonian. The time-independent electronic Schr¨ odinger equation describing the motion of electrons in the electrostatic field of the stationary nuclei comprises only term one, three and four of Eq. (2.54), where the electronic Hamiltonian and the electronic wave functions depend only parametrically on the coordinates of the nuclei. 2.3.1 DENSITY FUNCTIONAL THEORY The density functional theory (DFT) is a remarkable theory that allows one to replace the complicated N-electron wave function by a much simpler electron density ρ(~r). The basis for DFT is the proof by Hohenberg and Kohn [81] that the ground-state electronic energy is determined completely by the electron density ρ(~r)[79, 80, 82, 83]. The Hohenberg-Kohn says that the molecular energy, the wave function and all other molecular properties can be computed from the ground state electron probability density. Even with the proof that every specific electron density gives a specific energy (for the ground state), the connection between the energy and the functional has not been discovered. The hypothetical functional can be divided into three parts: the kinetic energy T[ρ(~r)], the core-electron attraction Ene[ρ(~r)], and the electron-electron repulsion Eee[ρ(~r)]. The electron-electron repulsion is usually divided into an exchange part K[ρ(~r)] and a Coulomb part J[ρ(~r)]. The exact kinetic energy (T) for a ground state is given by T= ∞ X i nihΨi| −1 2∇2 i|Ψii(2.55) 2.3. Quantum Chemistry 53 where Ψiare natural spin orbitals and nitheir occupation number, 0 ≤ni≥1. For an interacting system, Eq. (2.55) contains an infinite number of terms, making an exact solution unattainable. Kohn and Sham [81] presented a formalism for treating the kinetic part of the energy. In Hartree’s model [79, 80], the electrons move in an effective potential υH(~r)created by the atomic cores and a mean field created by the other electrons is given as: υH(~r) = − nuclei X a Za |~ Ra−~r |+Zρ(~r0) |~r −~r0|d~r0(2.56) where aruns over all the nuclei, Zais the nuclear charge on atom a,~ Radenotes the position of the nuclei and ρ(~r)is the electron density. In this approximation, a one-particle Schr¨ odinger equation can be defined as: −1 2∇2 i+υH(~r)Ψi=EiΨi(2.57) where Eiis the energy of the system. The mean density is given by ρ(~r) = N X i |Ψi(~r)|2(2.58) For the ground state, a summation over the Nspin orbitals with lowest eigenvalues is performed. Eq. (2.56-2.58) are solved self-consistently. From a trial guess of the electron density, a potential is obtained by Eq. (2.56). This potential is used in Eq. (2.57). The one-electron orbitals obtained are put in to Eq. (2.58). The procedure is repeated until a set of converged one-electron orbitals are obtained. For a system of independent non-interacting electrons, the common one-body potential, υKS is assumed. −1 2∇2 i+υKS(~r)Ψi=EiΨi(2.59) The local one-body potential υKS gives, by definition, a non-interacting electron density that has the same form as Eq. (2.57). On comparing Eq. (2.55) and Eq. (2.59), the expression for the kinetic energy of the noninteracting electrons (Ts[ρ]), becomes 54 Computational methods Ts[ρ] = N X i=1 hΨi| −1 2∇2 i|Ψii(2.60) where the sum iagain runs over all orbitals. The lower index sin Ts[ρ], denotes the singleelectron equations. As the electrons do interact in reality, the expression is not exact, and the following inequality holds where T[ρ]is the exact kinetic energy Ts[ρ]≤T[ρ](2.61) The remaining part of the kinetic energy defines the correlation contribution: Tc[ρ] = T[ρ]−Ts[ρ]≥0(2.62) where Tcis included in a exchange/correlation term Exc. The Kohn-Sham equations can be solved analogously to the Hartree equations, with the difference that the potential in Eq. (2.57), υH, is replaced by υeff : υeff (~r) = υ(~r) + Zρ(~r0) |~r −~r0|d~r0+υxc(~r)(2.63) where υxc(~r)denotes the local exchange/correlation potential and υ(~r)corresponds to the external potential. Within the KS formalism, the energy of the ground state can be obtained from EDF T [ρ] = X i Ei+Exc [ρ]−Zυxc(~r)ρ(~r)dv −1 2Zρ(~r)ρ(~r0) |~r −~r0|(2.64) or more generally, divided into its components: EDF T [ρ(~r)] = Ts[ρ(~r)] + Ene [ρ(~r)] + J[ρ(~r)] + Exc [ρ(~r)] (2.65) An exact energy expression has been obtained from Eq. (2.65). The first term is the kinetic energy for the non-interacting electrons and the second term is the core-electron attraction and the third term represents the Coulombic repulsion functional, the electrostatic electronelectron repulsion and the final term is the exchange/correlation term. Except the last term exchange/correlation energy, all other terms are solved exactly. The new challenge was to find a solution to Exc. 2.3. Quantum Chemistry 55 DIFFERENT DFT MODELS Most DFT models differ in their expression for the exchange/correlation energy Exc. It is always split into separate exchange and correlation terms, Exand Ec. The exchange energy is “by definition” given as a sum of contributions from αand βspin densities, as exchange energy only involves electrons of the same spin. The kinetic energy, the nuclear-electron attraction and Coulomb terms are trivially separable. The correlation energy contains contributions from the interactions between all electrons: Ex(ρ) = Eα x(ρα) + Eβ x(ρβ)(2.66) Ec(ρ) = Eαα c(ρα) + Eββ c(ρβ) + Eαβ c(ρα, ρβ)(2.67) Where ραand ρβare densities of the αand βspin. The total density is the sum of the αand β contributions, ρ=ρα+ρβthe largest contribution to Exc comes from the exchange part Ex. THE LOCAL DENSITY APPROXIMATION In the Local Density Approximation (LDA) it is assumed that the density can be treated locally as a uniform electron gas or the electron density is assumed to be slowly varying in space ELDA xc [ρ] = Zρ(~r)εunif xc (ρ)(2.68) εunif xc gives the exchange/correlation energy per electron in a uniform electron gas. The analytical expression for the exchange energy can be obtained from the Dirac formula ELDA xc [ρ] = −CxZρ4/3(~r)d~r (2.69) where Cxis given as: Cx=3 43 π1/3 In the more general case, where αand βdensities are not equal, LDA has been replaced by the Local Spin Density Approximation (LSDA). ELSDA xc [ρ] = −21/3CxZρ4/3 α(~r) + ρ4/3 β(~r)d~r (2.70) 56 Computational methods For the LDA correlation energy, no exact solution is known but it has been determined by Monte Carlo methods for a number of different densities. In order to use these results in DFT calculations, a suitable analytic interpolation formula is desirable. Vosko, Wilk and Nusari [84] fitted various data to obtain reliable analytic expressions for the LDA correlation energy and in general considered to be a very accurate fit. This exchange functional is used in the calculations. THE GENERALIZED GRADIENT APPROXIMATION The Generalized Gradient Approximation (GGA) [79, 80] is the improvement over the LSAD approach where the electron distribution is considered as a non-uniform electron gas. The idea behind this is not only considering the electron density, but also the gradient of the density, that is the change in the density at a given point. The term GGA was probably first introduced in connection with the PW86 functional, developed by Perdew and Wang in 1986 [85]. Perdew and Wang proposed modifying the LSDA exchange expression to that shown in Eq. (2.71) EP W 86 x=εLDA x(1 + ax2+bx4+cx6)1/15 (2.71) x=| ∇ρ| ρ4/3(2.72) where xis a dimensionless gradient variable, and a,band cbeing suitable constant where the summation over equivalent expression for αand βdensities is implicitly assumed. A general GGA functional has the form: εGGA =εGGA(ρα, ρβ,∇ρα,∇ρβ)(2.73) Several GGA functionals, for both exchange and correlation, have been proposed in the literature [85–90]. Becke [87] proposed a widely used correlation (B or B88) to the LSDA exchange energy: εB88 x=εLDA x+ ∆εB88 x(2.74) ∆εB88 x=−βρ1/3x2 1+6βx sinh−1x(2.75) The βparameter is determined by fitting to known atomic data and xis defined in Eq. (2.72). 2.3. Quantum Chemistry 57 Perdew and Wang have proposed an exchange functional similar to B88 to be used in connection with the PW91 correlation functional given below. εP W 91 x=εLDA x 1 + xa1sinh−1(xa2)+(a3+a4e−bx2)x2 1 + xa1sinh−1(xa2) + a5x2!(2.76) where aiand b are suitable constants and xis defined in Eq. (2.72). where xis given in Eq. (2.72) This GGA functional is used for the calculations. HYBRID METHODS The Kohn-Sham extensions revealed that improved expressions for electron exchange and correlation can lead to better wave functions and electronic energy. Therefore major attempts were carried out in the past to find improved and better electron exchange and correlation expressions. Wave function theory (Hartree-Fock) can in principle provide the exact exchange energy. EHF x=−1 2 n X i n X jZ Z Ψi(~r1)Ψj(~r1)Ψi(~r2)Ψj(~r2) |~r1−~r2|d~r1~r2(2.77) It could be assumed that the ideal exchange/correlation energy would then be obtained from Exc =EHF x+EDF T c(2.78) In hybrid DFT methods in general, the exchange term of DFT is corrected by a contribution from exact exchange energy. In the above, the DFT exchange has been replaced by HF exchange. BECKE3 The most popular hybrid methods are based on Becke’s three functional B3 [91]. It contains three adjustable parameters a0,axand ac: EB3 xc =a0EHF x+ (1 −a0)ELSDA x+ax∆EB88 x+ELSDA c+ac∆EGGA c(2.79) Becke used PW91 for the correlation part [86, 88, 92]. The name of the functional , B3PW91, implies its use of a three-parameter scheme, as well as GGA exchange and correlation functionals B and PW91, respectively. EB3LY P xc =a0EHF x+ (1 −a0)ELSDA x+ax∆EB88 x+ELSDA c+ac∆EP W 91 c(2.80) 64 Charge distribution in peptide conformations effect of these methods on dipole moments. They checked the ability of these charges by reproducing the experimental dipole moments for 89 molecules. Various methods like RESP and CHELP-Bow were used to drive partial atomic charges of the disaccharides from the quantum chemical electrostatic potential and moments [137]. The atomic charges for various conformations was studied by Richards and co-workers [138] who reported two related methods, both of them employ a constrained minimization of the difference between the quantum mechanical and classical MEP (Molecular Electrostatic Potential) with respect to the atomic charges. The first method involves determining the MEP and constraining the charges to reproduce the dipole at an alternative geometry. The second method involves determining the MEP for each conformation of interest and and weighting the MEP for each conformation according to the appropriate Boltzmann factor. The charge distribution of the tripeptide Ala-Asn-Ala with different conformations have been analyzed using the hybrid density functional method (B3LYP) with 6-31G* basis set. The population analysis methods such as Mulliken, NPA and the electrostatic potential methods CHELPG, MK and RESP are used to analyze the charge distribution of the tripeptide Ala-AsnAla. The extensive studies on the charge methods can provide a better understanding of each method and the parameters which influence the partial atomic charges. The different charge methods discussed in this chapter are used to obtain the point charges for the electrostatic calculations and for the calculation of pKavalues of CuBligands in aqueous solution. 3.1 THEORETICAL BACKGROUND In this section the theory behind charge methods, the advantage and disadvantages of each charge methods are explained. 3.1.1 WAVE FUNCTION BASED METHODS MULLIKEN POPULATION ANALYSIS In Mulliken population analysis, the electrons are divided up amongst the atoms according to the degree to which different atomic orbital (AO) basis functions contribute to the overall wave function [79, 80, 139, 140]. The total number of electrons are obtained by expanding the wave function in its AO basis set, N= electrons X jZψj(~rj)ψj(~rj)d~rj(3.1) = electrons X jX r,s Zcjrϕr(~rj)cjsϕs(~rj)d~rj(3.2) 3.1. Theoretical background 65 = electrons X jX r c2 jr +X r6=s cjrcjsSrs(3.3) where ψis the wave function of molecule j,rand sare the index of AO basis function ϕ,cjr is the coefficient of the basis function rin molecular orbital j,cjs is the coefficient of the basis function sin molecular orbital j.Srs is the usual overlap matrix element [79]. From Eq. (3.3) its clear that the total number of electrons are divided into two sums, one including only the squares of single AO basis functions and the other including products of two different AO basis functions. The first term in the parentheses in Eq. (3.3), is the sum of electrons associated with single basis function only and the second term, which represents the electrons ‘shared’ between basis functions. The total number of electrons per atom Nkis given by Nk= electrons X jX r∈k c2 jr +X r,s∈k,r6=s cjrcjsSrs +X r∈k,s/∈k cjrcjsSrs(3.4) The orthonormality of basis functions of different angular momentum both residing on the same atom kcauses many terms in the second sum of Eq. (3.4) to be zero. The Mulliken partial atomic charge is defines as qk=Zk−Nk(3.5) where Zkis the nuclear charge of atom kand Nkis computed according to the Eq. (3.4). With minimal or small split-valence basis sets, the Mulliken charges tend to be reasonably intuitive. The use of non-orthogonal basis set in the Mulliken analysis, however, can lead to some undesirable results. In Mulliken analysis where the total number of electrons are divided over AO basis functions, it is possible for individual basis functions to have occupation numbers greater than 1 or less than 0. Such situations obviously can have no physical meaning. The electrons are also divided equally in the case atoms with different electronegativities. Mulliken partial charges prove to be very sensitive to basis-set size, so that comparisons of partial atomic charges from different levels of theory are in no way possible. Moreover, with very complete basis sets, the Mulliken charges have a tendency to become unphysical. NATURAL POPULATION ANALYSIS Natural Atomic Orbitals (NAOs) were defined by P. O. L¨ owdin [141] as the eigenfunctions of the first-order reduced density matrix [142]. They are intrinsic orbitals {θk}that arise as eigenfunctions of the first-order reduced density operator ˆ Γ, ˆ Γθk=qkθk(3.6) 66 Charge distribution in peptide conformations which is formed by ‘reducing’ the wavefunction probability distribution to the single-particle level, ˆ Γ = NZψ(1,2, ..., N)ψ∗(10,2...N)dτ2...dτN(3.7) and whose eigenorbitals are hence ‘natural’ to the ψitself. The subsystems of the N-particle system are best formulated in terms of reduced density operators. In particular, the squared amplitude |hψ(1,2,...,N)|φ(1)i|2that an electron of ψ(1,2,...,N) is ‘in’ orbital φ(i.e., the population of φ(1) in the wavefunction) is rigorously expressed, for any possible orbital φ, as qφ=hφ|ˆ Γ|φi(3.8) The occupancies qφare intrinsically non-negative and limited by the Pauli exclusion principle e.g., for the spatial orbital φ(r), 0≤qφ≤2(3.9) The sum of occupancies qkover any complete orthonormal set {θk}accounts for all Nelectrons, X k hφk|ˆ Γ|φki=X k qk=N(3.10) Mulliken population commonly fail to satisfy the physical constrains (Eq. 3.9 and 3.10). The natural orbitals {θk}of Eq. (3.6) are intrinsically best possible for considering maximum occupancy in the smallest possible number of orbitals. The starting point for the Natural Population Analysis (NPA) is the determination of an optimal set of effective atomic orbitals for the molecular environment. To find the maximum-occupancy orbitals for atom A, the local eigenvectors ˜ θ(A) iof the onecenter atomic block ˆ Γ(A) associated with basis AOs on atom A is obtained first. ˆ Γ(A) ˜ θ(A) i= ˜q(A) i˜ θ(A) i(3.11) After the maximum-occupancy orbitals {˜ θ(A) i}have been obtained for each atom A, it is necessary to restore full orthogonality. The ‘pre-orthogonal’ orbitals {˜ θ(A) i}are then subjected to an occupancy-weighted symmetric orthogonalization to find final NAOs {θi(A)}. Once the final orthonormal NAOs θi(A) are obtained, the Eq. (3.8) are used to rigorously compute their ‘natural’ populations qi(A), 3.1. Theoretical background 67 qi(A) =hθ(A) i|ˆ Γ|θ(A) ii(3.12) which may be summed to give the total number of electrons NA NA=X i qi(A) (3.13) and ‘natural charge’ QAon atom A with atomic number ZAis given as QA=ZA−NA(3.14) In accord with Eq. (3.9 and 3.10), the natural populations qi(A) necessarily satisfy the Pauli principle (0 ≤qi(A) ≤2) and sum correctly to the total number of electrons. The Eq. (3.12-3.14) summarize ‘natural population analysis’. The “Natural Population Analysis” is an alternative to the conventional Mulliken population analysis [79, 80, 143, 144]. It seems to exhibit improved numerical stability and better in describing the electron distribution in compounds of high ionic character. The most appealing feature of the NPA scheme is that each atomic partial charge effectively converges to a stable value with increasing basis-set size. 3.1.2 POTENTIAL BASED METHODS Widely distributed methods to generate atomic partial charges from QC computations are called potential-based methods such as CHELP, CHELPG, MK and RESP method. The charges are calculated by matching the electrostatic potential (ESP) from the instantaneous atomic partial charge distribution with the QC ESP computed on optimized geometries in vacuum. These methods have in common that they use values of the electrostatic potential on spatial positions that are distributed around the atoms of the molecule. All potential-based methods take into account only the potential points outside the vdW radii of the atoms to determine charges. As they cannot be determined undoubtedly, effective vdW radii are a matter of debate in the literature and are treated differently in the potential-based methods. Different schemes are currently in use to select the points outside the vdW radii of the atoms of the molecule. MOLECULAR ELECTROSTATIC POTENTIAL A truncated multipole expansion is a computationally convenient single-center formalism that allows one to quantitatively compute the degree to which a positive or negative test charge is attracted or repelled by the molecule that is being represented by the multipole expansion [79]. The MEP can be computed exactly on any given position ~r 68 Charge distribution in peptide conformations VMEP (~r) = nuclei X k Zk |~r −~rk|−ZΨ(~r0)1 |~r −~r0|Ψ(~r0)d~r0(3.15) Where Zkis the nuclear charge and the sum run overs all the atom k. The MEP is an observable, although in practice it is rather difficult to design appropriate experiments to measure it. Computationally, the molecular electrostatic potential is usually evaluated as VMEP (~r) = nuclei X k Zk |~r −~rk|−Zρ(~r0) |~r −~r0 |d~r0(3.16) where ρis the electron density. The partial atomic charges are typically derived from MEP. The common notation MEP is replaced with ESP for, ‘electrostatic potential’. The ESP is the most obvious property to obtain the partial atomic charges that will be useful to model molecular-molecular interactions at short and long range. All ESP charge-fitting schemes involve determining atomic partial charges qk, VESP (~r) = nuclei X k qk |~r −~rk|(3.17) Although there are many positive aspects of the ESP derived charges, often they are found to be very sensitive to small changes to the molecular structure or in the sampling of the ESP. CHELP METHOD In CHELP (CHarges from ELectrostatic Potentials) method [79, 145] the electrostatic potential is determined at a selected number of points around the molecule chosen in spherical shells, 1˚ A apart, of fourteen nearly symmetrically placed points around each atom. Point which fall within the vdW radius of any of the atoms are discarded due to the large distortions caused by close proximity to the nuclei. Typically 100-300 points are selected in the region extending to 3 ˚ A from the vdWs surface of the molecule is considered. From Eq. (3.16) the ESP is calculated quantum chemically and a least square procedure is used to fit the charge qjto each atomic center jin the molecule. The calculated classical ESP is given by Eq. (3.17). The χ2 esp to be minimized by least-squares procedure is defined as in Eq. (3.18). χ2 esp =X i (Vi−Ei)2(3.18) At the minimum, for all j 3.1. Theoretical background 69 ∂(χ2 esp) ∂qj = 0 (3.19) where ∂(χ2 esp) ∂qj =−2X i Vi−Ei rij = 0 (3.20) and thus a system of equations can be formed and solved in matrix form. The best least square fit of charges to the electrostatic potential is obtained by finding the minimum of y, y(q1, q2, ...qn) = m X i=1 [Vi−Ei(q1, q2, ...qn)]2(3.21) where mis the total number of points to fit, Viis from QC wave functions as given in Eq. (3.16) and Eiis the electrostatic potential in the monopole approximation and given in Eq. (3.17). The minimum of ycan be obtained by finding the stationary points of the Lagrangian function z: z(q1, q2, ...qn) = y(q1, q2, ...qn) + λg(q1, q2, ...qn)(3.22) where gis the constraint imposed on the fit. For example, following constraints can be used: constraining the total charge, molecular dipole, the charges of a group of atoms to a particular value and tow or more atoms having equivalent charges. The linear least square problem with imposed constraints can be solved using the Lagrangian multiplier method. The function corresponding to total charge constraint is: g(q1, q2, ...qn) = Xqi−qtot = 0 (3.23) The extremum points of zare then obtained by solving ∂z/∂λ = 0 and ∂z/∂qk= 0. A set of n+1 equations are obtained where nis the number of atoms. The solution of corresponding matrix equation yields the charges. The CHELP method has rotational invariance problems due to the sparse point sampling and dependency upon the molecular coordinate system for choosing test points. CHELPG METHOD The principle difference between the CHELPG (CHarges from ELectrostatic Potentials using a Grid based method) and CHELP is that CHELPG [79, 146] procedure employs a point selection 70 Charge distribution in peptide conformations algorithm based upon regularly spaced points. The algorithm involves defining a cube of points spaced 0.3-0.8 ˚ A apart containing the molecule and including 2.8 ˚ A of head space on all sides. Figure 3.1 describes the situation for water as a simple example. All points outside the 2.8 ˚ A radius are also eliminated from the fitting procedure. The remaining points form a relatively homogeneous layer around the molecule. The electrostatic potential at each sample points is calculated analytically from the wave function and geometry data. These data are then used as input for the Lagrange multiplier least-squares routine, which has been constrained to fit the exact molecular charge, but not to the molecular dipole moment. CHELPG typically utilize 5000-6000 points. The CHELPG charges are obtained for all the five conformations by varying the grid points per square angstrom. In the present study the grid points in the CHELPG method is varied from 1 to 6 (points per square angstrom) in ESP fit. O HH 28 pm 28 pm Figure 3.1. CHarges from ELectrostatic Potentials using a Grid based method. The dimensions of the cube are chosen such that the molecule is located at the center of the cube, adding 28.0 pm head space between the molecule and the end of the box in all three dimensions. One of the weak points of CHELPG is the treatment of larger systems, in which some of the innermost atoms are located far away from the points at which the ESP is evaluated. In such situation, variations of the innermost atomic charges will not lead to significant changes in the ESP outside the molecule and fitting of these atomic charges will therefore not result in meaningful results. Representative atomic charges for flexible molecules should therefore be derived as average values over several conformers. MK METHOD The Merz Kollman (MK) method [79, 147, 148] uses the Connolly surface algorithm [149] to calculate a shell or a number of shells of points at a specified distance from the molecule and then calculate the electrostatic potential at these points. The use of successive shells of various multiples of the vdWs radii, ranging from 1.4 to 2.0 times the radius (smaller radii lead to positive potentials in the lone-pair region). The shells of 1.4, 1.6, 1.8, and 2.0 3.1. Theoretical background 71 times the vdWs radius are used to avoid artifacts from being too close to the atoms. The density of points in these surfaces depended on the size of the molecule and typically 200300 points/molecules are generated by varying 1-5 points per square angstrom. The Merz Kollman method is described for water in Figure 3.2. The water molecules is covered with four shells with 1.4, 1.6, 1.8, 2.0 and 2.2 times the vdWs radius. In the present study the MK charges are obtained for all the five conformations by varying the grid points and the shells. The grid points varied from 1 to 6 and the shell is varied from 2 to 9 to the maximum limit with 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8 and 3.0 times the vdWs radius. Vdw envelope HH O Shell 1 (scaling 1.4) Shell 2 (scaling 1.6) Shell 3 (scaling 1.8) Shell 4 (scaling 2.0) Shell 5 (scaling 2.2) Figure 3.2. Merz-Singh-Kollman method. The ESP is calculated at a number of grid points located on several shells around the molecule. The shells are constructed as an overlay of van der Waals spheres around each atom. The shells of 1.4, 1.6, 1.8,2.0 and 2.2 times the vdWs radius are shown. The major strength of the electrostatic potential charges is that they optimally reproduce the intermolecular interaction properties of the molecules with a simple two-body additive potential, if suitably accurate level of quantum chemical calculation is used to derive the ESP around the molecule. The major weakness of these charges is they are not easily transferable between common functional groups in related molecules [133]. RESTRAINED ELECTROSTATIC POTENTIAL METHOD (RESP) In RESP method [150], the restraints are introduced in the form of penalty function into the fitting procedure considerably to overcome the limitations of other ESP charge methods. With the charge fitting procedure χ2 esp (see section 3.1.2), an addition penalty function χrstr is added. So χ2is to be minimized now χ2=χ2 esp +χ2 rstr (3.24) 72 Charge distribution in peptide conformations and the least squares minimum is defined as ∂(χ2) ∂qn =∂(χ2 esp) ∂qn +∂(χ2 rstr) ∂qn (3.25) The simple harmonic Eq. (3.26) or hyperbolic Eq. (3.27) function is used as a penalty function. χ2 rstr =aX m (q0n−qn)2(3.26) χ2 rstr =aX m ((q2 n+b2)1/2−b)(3.27) where ais the scaling factor determining the strength of the restraint and q0is the target for the restraint. Mulliken charges are used as a target charges [150]. Generally the RESP charges are used in molecular mechanics form any organic and bio-organic systems. In the present study, the RESP charges are obtained by using both Mulliken and NPA charges as initial charges for the harmonic restraint. The RESP charges are obtained from both electrostatic potentials calculated from CHELPG and MK methods. MULTIPOLE DERIVED CHARGES Another set of charges used in the work is multipole derived charges. The multipole derived charges (MDC) are used to obtain the solvation energies based on FDPB method combination with DFT level of theory (see section 5.3.1). The MDC analysis is based on the ideas used for the Dipole Preserving Charge analysis [151] but formulated by using more accurate atomic dipoles. There are three stages involved in the method: First the molecular electron density is written as a sum of atomic densities. Second, from these atomic densities a set of atomic multipoles are defined, which is used to get the electrostatic potential outside the charge distribution. Third, these atomic multipoles are reconstructed by using a scheme that distributes charges over all atoms to reproduce these multipoles exactly. It is claimed that this method can overcome the drawbacks of the ESP methods because this method uses the atomic multipole expansion over molecular multipole expansion. A detailed description of atomic multipole expansion method is given in Ref [152]. The MDC used in the present work are charges derived from monopole (MDCm), dipole (MDCd) and quadrupole (MDCq) moments. 3.2 STRUCTURE PREPARATION AND COMPUTATIONAL METHODS The initial structure of the tripeptide Ala-Asn-Ala was builded using CHARMM program [153] package. The C-terminal alanine and the N-terminal alanine are capped with acetelate and methylate groups respectively. Five different conformations are generated by changing the 3.2. Structure preparation and Computational Methods 73 dihedral angles of atoms 6, 7, 17 and 18 from 120◦to 180◦by step of 15◦. The tripeptide considered for the charge calculation is shown in Figure 3.3 1Ala − 2Asn − 3Ala N CA O N CA CG OD1 C N CA NT C CAT CB ND2 1 2 3 4 5 6 7 8 9 11 12 13 15 10 16 17 18 19 21 22 23 CAY CY CB C OY CB O O 14 20 Figure 3.3. The tripeptide Ala-Asn-Ala considered for the charge calculation is shown. The tripeptide Ala-Asn-Ala is shown with the PDB atom type. Each conformation was optimized by fixing the dihedral angle. All the conformations were optimized by hybrid density functional theory B3LYP with 6-31G* basis set. The geometry optimizations were performed using GAUSSIAN 03 program [154]. The CHELPG and MK charges were also obtained using the same program package. The RESP charges were obtained using RESP program [155]. The complete set of calculations performed on tripeptide Ala-Asn-Ala are shown in Figure 3.4. Population analysis were performed for all the conformations. The choice of the population analysis method includes the traditional Mulliken population analysis [139, 140] and NPA [143, 144], both of these methods are based on the molecular wave function, but NPA uses the natural atomic basis function to overcome the drawback of the Mulliken population analysis. Various methods are available to fit the molecular electrostatic potential to the grid points with fitting procedures. The three methods considered in our studies are the CHELPG [145] , MK [147, 148], and RESP [156] methods. The CHELPG method is a grid based method where points are selected in a regularly spaced cubic grid. To obtain charges independent of the conformations of the tripeptide, the sampling points were increased by varying the grid points in CHELPG and MK method and the charges were calculated for each conformation. In CHELPG method, the grid points per square angstrom are varied from 1 to 6 in the ESP fit by invoking the GAUSSIAN 03 option IOp(6/42=N) where N is the points per square angstrom. MK method uses the Connolly surface algorithm [149] to calculate a shell or a number of shells 80 Charge distribution in peptide conformations MODEL 3 In model 3, the ESP obtained from MK method by varying grid points from 1-6 are used to calculate RESP charges and the Mulliken charges are used as initial charges for harmonic restraint. As in model 1, when less grid points and strong restraint 0.5 or 0.1 are used in the first stage fit, the RESP charges approach the Mulliken charges. When the grid points are increased from 3 to 4, there in no change in the charge and increasing the sampling points do not have any effect on the charges. The calculations show that at least 3 point per square angstrom should be used to get a stable charges and the restraint factor of 0.0005 should be used to get a reasonable charges. MODEL 4 In Model 4, the NPA charges are used as initial charges with the combination of ESP from MK method. The trends are almost same. Analysis of the RESP charges also show that more than 3 points per square angstrom should be used to get stable charge with combination of restraint factor 0.0005. MODEL 5 In Model 5, the RESP charges are calculated by using the ESP obtained from MK method by varying the shells from 2 to 9. There is an irregular trend in the RESP charges with respect to the conformations. The analysis shows that at least more than 4 shells should be used to get a stable charges for each conformations. All the models show that the RESP charges approach the initial charges (Mulliken and NPA) when strong restraint is used and the RESP charges approach the ESP charges (CHELPG and MK) when weak restraint is used. The restraint beyond 0.00005 does not make any change in the charges. 3.4 CONCLUSIONS The Mulliken, NPA, CHELPG, MK and RESP atomic charges are derived for different conformation of the tripeptide Ala-Asn-Ala. All the partial atomic charges were obtained from hybrid density function theory B3LYP with 6-31G* basis set. The sensitivity of charges with different conformations and the influence of wave-function based (Mulliken and NPA) and electrostatic potential based charge methods (CHELPG, MK and RESP) have been explored in detail. The factors influencing the charges are varied which includes increasing the points per square angstrom both in CHELPG and MK method and the shells are varied from 2-9 the maximum limit in MK method. The RESP charges were calculated by using the electrostatic potential from CHELPG and MK method using Mulliken and NPA charges as initial charges with harmonic function as penalty function. The RESP charges were also obtained by using different scaling factor from a strong to weak restraint. The atomic charges obtained from both CHELPG and MK method are dependent on the conformations. The Mulliken and NPA charges show less conformational dependency compare to 3.4. Conclusions 81 the CHELPG and MK charges. Difference in points selection for the ESP fits can dramatically affect the partial charges. The magnitudes of the calculated point charges were found to be somewhat dependent on the choice of points per square angstrom and the shells. The magnitude of the charges are more or less same for the both CHELPG and MK method when grid points of 3 to 4 per square angstrom were used. The results show that minimum 2 to 3 grid points per square angstrom should be used in both methods and increasing the grid points to maximum limit 6 does not change the magnitude of charge much. The MK charges show less dependent on conformation when minimum 3 to 4 shells were used and the charges remains almost same when shells are increased to maximum limit 9. The results show that minimum 3 to 4 shells should be used in MK method. The RESP charges are mainly dependent on the strong or weak restraint used. When strong restraints are used the charges are closer to the Mulliken or NPA charges and if weak restraints are used the charges are close to the ESP charges. The magnitude of the charges exhibit much smaller fluctuations when restraint factor of 0.0005 was used in the first stage fit. The restraint weight of 0.0005 appears to be optimal for obtaining RESP charges for tripeptide Ala-Asn-Ala. The model 3 combination works fine to obtain RESP charges for tripeptide Ala-Asn-Ala. The charges predicted by the Mulliken population analysis and NPA schemes are larger than those predicted by the ESP schemes. This conclusions should be carefully reconsidered for other systems with polar bonds in which electrons correlation effects may be more significant. 82 Charge distribution in peptide conformations A.1 A.2 A.3 B.1 B.2 B.3 Figure 3.6. Atomic charges for heavy atoms in the backbone and side chains for different conformations of tripeptide Ala-Asn-Ala obtained from CHELPG and MK methods are shown in the plots. (A.1). CHELPG charges for heavy atoms in the backbone and side chain where 1 point per square angstrom is used in the CHELPG method. (A.2). CHELPG charges for heavy atoms in the backbone and side chain where 2 points per square angstrom are used in the CHELPG method. (A.3). CHELPG charges for heavy atoms in the backbone and side chain where 3 points per square angstrom are used in the CHELPG method. (B.1). MK charges for heavy atoms in the backbone and side chain where 1 point per square angstrom is used in the MK method. (B.2). MK charges for heavy atoms in the backbone and side chain where 2 points per square angstrom are used in the MK method. (B.3). MK charges for heavy atoms in the backbone and side chain where 3 points per square angstrom are used in the MK method. The black circle (◦) corresponds to conformation 1 (1-20), red square (2) corresponds to conformation 2 (2-135), green diamond () corresponds to conformation 3 (3-150), blue triangle up (M) corresponds to conformation 4 (4-165), magenta star (∗) corresponds to conformation 5 (5-180) 3.4. Conclusions 83 A.4 A.5 A.6 B.4 B.5 B.6 Figure 3.7. Atomic charges for heavy atoms in the backbone and side chains for different conformations of tripeptide Ala-Asn-Ala obtained from CHELPG and MK methods are shown in the plots. (A.4). CHELPG charges for heavy atoms in the backbone and side chains where 4 points per square angstrom are used in the CHELPG method. (A.5). CHELPG charges for heavy atoms in the backbone and side chains where 5 points per square angstrom are used in the CHELPG method. (A.6). CHELPG charges for heavy atoms in the backbone and side chains where 6 points per square angstrom are used in the CHELPG method. (B.4). MK charges for heavy atoms in the backbone and side chain where 4 points per square angstrom are used in the MK method. (B.5). MK charges for heavy atoms in the backbone and side chains where 5 points per square angstrom are used in the MK method. (B.6). MK charges for heavy atoms in the backbone and side chains where 6 points per square angstrom are used in the MK method. The black circle (◦) corresponds to conformation 1 (1-20), red square (2) corresponds to conformation 2 (2-135), green diamond () corresponds to conformation 3 (3-150), blue triangle up (M) corresponds to conformation 4 (4-165), magenta star (∗) corresponds to conformation 5 (5-180) 84 Charge distribution in peptide conformations C.1 C.2 C.3 C.4 C.5 Figure 3.8. Atomic charges for heavy atoms in the backbone and side chains for different conformations of tripeptide Ala-Asn-Ala obtained from CHELPG and MK methods is shown in the plots. (C.1). MK charges for heavy atoms in the backbone and side chains of conformation 1 (1-120) where shells varied from 2-9 and 4 points per square angstrom is used in the MK method. (C.2). MK charges for heavy atoms in the backbone and side chains of conformation 2 (1-135) where shells varied from 2-9 and 4 points per square angstrom is used in the MK method. (C.3). MK charges for heavy atoms in the backbone and side chains of conformation 3 (1-150) where shells varied from 2-9 and 4 points per square angstrom is used in the MK method. (C.4). MK charges for heavy atoms in the backbone and side chains of conformation 4 (1-165) where shells varied from 2-9 and 4 points per square angstrom is used in the MK method. (C.5). MK charges for heavy atoms in the backbone and side chains of conformation 5 (1-180) where shells varied from 2-9 and 4 points per square angstrom is used in the MK method. 3.4. Conclusions 85 D.1 D.2 D.3 D.4 D.5 Figure 3.9. Atomic charges for heavy atoms in the backbone and side chains for different conformations of tripeptide Ala-Asn-Ala obtained from RESP methods is shown in the plots. (D.1). RESP charges for heavy atoms in the backbone and side chains of conformation 1 (1-120). (D.2). RESP charges for heavy atoms in the backbone and side chains of conformation 2 (1-135). (D.3). RESP charges for heavy atoms in the backbone and side chains of conformation 3 (1-150). (D.4). RESP charges for heavy atoms in the backbone and side chains of conformation 3 (1-165). (D.5). RESP charges for heavy atoms in the backbone and side chains of conformation 3 (1-180). The black circle (◦) corresponds to RESP charges where ESP from CHELPG method is used with combination of Mulliken charges as initial charges. The red square (2) corresponds to to RESP charges where ESP from CHELPG method is used with combination of NPA charges as initial charges. The green diamond () corresponds to RESP charges where ESP charges MK method is used with combination of Mulliken charges as initial charges. The blue triangle up (M) corresponds to RESP charges where ESP from MK method with combination of NPA charges as initial charges. 86 Charge distribution in peptide conformations CHAPTER 4 PROTONATION AND REDOX POTENTIALS OF CYTOCHROME cNITRITE REDUCTASE Cytochrome coxidase is a large membrane protein designed to utilize the energy of electron transfer and oxygen reduction to pump protons across the membrane. The molecular mechanism of the energy conversion process is not been well understood [8, 40–42]. Electrostatic calculations on other proteins with simpler, better resolved structures can help to understand the possible mechanism of electron transfer in cytochrome coxidase. Cytochrome cnitrite reductase can serve as such a simple system. Cytochrome cnitrite reductase catalyzes the six electron reduction of nitrite to ammonia. This second part of the respiratory pathway of nitrate ammonification is the key step in the biological nitrogen cycle [158–161]. Cytochrome cnitrite reductase is a multiheme cenzyme and its active site is a protoporphyrin IX which is covalently linked to the protein backbone. A lysine was found to replace the usual histidine as a proximal ligand to the heme iron [162–164]. The five hemes in the monomer of cytochrome cnitrite reductase are in close contact with Fe-Fe distance of 9 to 12.8 ˚ A. The cytochrome cnitrite reductase dimer is shown in Figure 4.1 and the arrangement of hemes in the dimer is shown in Figure 4.2. The Heme I forms the active site and the other hemes, Hemes II, III and IV are almost coplanar with the catalytic Heme I. Hemes II and V are farther apart and are not coplanar with Hemes I, III and IV. All the hemes except Heme I are bis-histidinylcoordinated and linked to the peptide backbond by thioether bonds to the cysteine residues of a classical heme-binding motif for periplasmic proteins, Cys-X1-X2-Cys-His. Heme I however has the binding motif Cys-X1-X2-Cys-Lys, where the nitrogen atom of the lysine replacing that of histidine. The cytochrome cnitrite reductase catalyzes the reduction of both nitrite and sulfide with high specific activity and gets electrons from the membranous quinol pool, thereby generating a proton motive force. Heme I is clearly the substrate-binding site and previous spectroscopic studies [165, 166] in a similar nitrite reductase had proved that one of the heme was highspin and it was the active site of the enzyme. The cytochrome cnitrite reductase can be found in the periplasm and forms a stable, membrane associated complex with its electron donor NrfH, a member of the NapC/NirT family of the tetraheme cytochromes [167, 168]. Cytochrome cnitrite reductase possesses more than 180 protonatable groups, of which ten of them are the heme propionates in each monomer. In total, cytochrome cnitrite reductase 87 88 Protonation and redox potentials of cytochrome cnitrite reductase contains 10 c-type hemes. In ten hemes, eight hemes are c-type heme i.e., two histidines coordinating to the iron and in other hemes the lysine is covalently linked as a proximal ligand to the heme iron of the c-type heme. Although all c-type hemes are of the same chemical nature, their midpoint potentials differ. The cytochrome cnitrite reductase is a homodimeric enzyme with 10 c-type hemes which are arranged such that the nearest neighbors are in close proximity. It has been reported [164, 169] that the reduction potential of the 10 heme centers ranges from ca. -30 to -320 mV in cytochrome cnitrite reductase from Escherichia coli. The protein film voltammetric experiments revealed that the heme oxidation state has a profound, and often unanticipated effect on the interactions with substrate molecules, nitrite, hydroxyl amine and inhibitor cyanide [169]. The oxidation probabilities of the hemes and protonation probabilities of the heme propionates and redox potentials of the hemes in cytochrome cnitrite reductase are calculated to study the redox potential profiles. The electrostatic interactions between the protonatable and redox-groups are studied by solving the Poisson-Boltzmann equation using finite difference method. The electrostatic calculations are performed in cytochrome cnitrite reductase, a less complex protein before performing calculations on the complex protein like cytochrome c oxidase which contains four redox centers CuA, heme a, heme a3and CuB. The aim of the present study is to obtain the redox mid-point potentials of the hemes to understand the redox potential profiles of hemes when the catalytic Heme I was blocked with the inhibitor cyanide. Figure 4.1. The nitrite reductase dimer and the heme arrangement. A front view with the dimer axis orientated vertically, the five hemes in each monomer (stick model) and the Ca2+ (yellow) are shown. 4.1. Hemes and the calcium binding site in cytochrome cnitrite reductase 89 Heme I Heme III Heme II Heme IV Heme V Figure 4.2. The arrangement of hemes in nitrite reductase dimer is shown. The overall orientation of hemes corresponds to nitrite reductase dimer is shown. Hemes in the left monomer are labeled according to their attachment to the protein chain. The same labels are used throughout the study. 4.1 HEMES AND THE CALCIUM BINDING SITE IN CYTOCHROME c NITRITE REDUCTASE Hemes I, III and IV are close enough to allow direct π-electron interaction of the porphyrin rings. The propionate side chains of heme I form part of the active-site cavity, while those of Heme IV are exposed to the solvent and the Heme III propionates are hydrogen-bonded inside the protein. All the hemes show slight distortion from ideal planarity and it is more pronounced in Heme II and less in Heme I. It is suggested [162] that Heme II could function as the entry point for electrons (see Figure 4.2). The Ca2+ binding site appears to be an essential structural feature in the overall architecture of the enzyme and the region surrounding the calcium binding site is one of the most highly conserved parts of the whole sequence [163]. The calcium binding site of cytochrome cnitrite reductase with the important residues are shown in Figure 4.3. The residues Tyr281, Lys274, Gln276 and His277 act as ligands for the calcium binding site. It was proposed that the calcium binding site hold the key residues needed for catalysis [162]. The set of tyrosine residues near the calcium binding site and active site might play a role in the reaction mechanism of cytochrome cnitrite reductase by forming possible radical intermediates of the stepwise reduction of nitrite to ammonia [170]. 96 Protonation and redox potentials of cytochrome cnitrite reductase Hemes Mid-point redox potentials in [mV] Oxd.aRed.bExp.c Heme II -210 -210 -323 (Heme IV) Heme III -30 -50 -37 (Heme II) Heme IV -270 -270 -107 (Heme III) Heme V -210 -210 -323 (Heme V) Table 4.3. The calculated mid-point redox potentials of Hemes II to V at pH=7. The mid-point redox potentials obtained for the Hemes II-V in the reduced and oxidized states of Heme I are given. The experimental mid-point redox potentials are arranged where the experimental values match closly with the calculated mid-point potentials. The order of assignment of redox-potential to a particular heme is very controversial [169]. aThe mid-point redox potentials were calculated for Hemes II-V by fixing the catalytic Heme I in the oxidised state. bThe mid-point redox potentials were calculated for Hemes II-V by fixing the catalytic Heme I in the reduced state. cThe experimental mid-point redox potentials of Hemes II-V from Ref [169]. Figure 4.6. Oxidation probabilities of different hemes in cytochrome cnitrite reductase from W. succinogenes at pH=7. The hemes do not show a standard Nernst titration behavior because of the interaction with each other and with the titratable groups in their vicinity. of these interactions, titration curves of titratable groups in proteins can deviate considerably from sigmoidal Nernst or Henderson-Hasselbalch titration curves of isolated titratable groups. Thus, it is not always possible to assign a unique mid-point potential or pKavalue to a specific titratable group. Therefore, the pH value at which the protonation probability of the protonatable group is 0.5 is often used to describe the titration behavior. Likewise, the 4.3. Results and Discussion 97 solution redox potential E◦ 1/2at which the oxidation probability of the redox-active group is 0.5 is used to describe the titration behavior. The oxidation probabilities of different hemes in cytochrome cnitrite reductase at pH=7 are given in Figure 4.6. The titration curves deviates from the standard Nernst titration curves. This deviation is expected because of the interaction between the hemes and the interaction of hemes with other protonatable groups and heme propionates. Heme II Heme III Heme IV Heme V Figure 4.7. Oxidation probabilities of the Hemes II-V. Calculated oxidation probabilities of heme II-V of cytochrome cnitrite reductase from W. succinogenes. The probabilities are color coded as indicated by the scale next to the diagram. The green color indicates that the oxidation probability of the hemes is about 0.5. The computed redox potentials for hemes range from –270 to –30 mV. However the mid-points cannot be compared with the experimental mid-point potential because the experimental midpoint potentials were obtained using nitrite as inhibitor. The redox potentials in the present studies were computed by blocking the catalytic Heme I by a cyanide where the blocking prevents the electron transfer from the catalytic heme to the inhibitor. Butt et al. [169] obtained redox potential of –107 mV for Heme III. The redox potential of –30 mV was obtained for Heme III from present calculations. On comparing the redox potentials of Hemes II-V, in 98 Protonation and redox potentials of cytochrome cnitrite reductase the oxidized and in the reduced state of Heme I, the redox potentials remain almost same expect for Heme III where the difference is –20 mV. Hemes II and V have same mid-point potential even though the distance between the hemes are more compare to heme IV and V. The Heme III has high redox potential compared to other hemes. The oxidation probabilities of the Hemes II-V are shown in Figure 4.7. Within the protein, the heme mid-point redox potentials are affected by charges on the ionizable amino acids, polar groups and other hemes. The mid-point potential of the Hemes II-V decreases with the increase in pH. The oxidation probabilities of the Hemes II-V (see Figure 4.7) show that redox states of the hemes are strongly depend on the pH. The pH dependence is less pronounced in Hemes II and V which are not coplanar with other hemes and more pronounced in Heme III and IV. PROTONATION PROBABILITIES OF HEMES I-V PROPIONATES a.1. b.1. Figure 4.8. Protonation probabilities of the hemes propionates A and D of Heme I. (a.1)-(b.1). Calculated protonation probabilities of heme propionates A and D of Heme I. The probabilities are color coded as indicated by the scale next to the diagram. The green color indicates that the protonation probability of the residues is about 0.5. The protonation probabilities of heme propionates (A and D) from Heme I and Heme II-V are shown in Figure 4.8 and 4.10 respectively. The interaction energies between the hemes and heme propionates are given in Figure 4.9. The protonation probabilities of Heme I propionates are strongly depend the redox states of other hemes. The protonation probabilities of Heme I propionate A (plot a.1) strongly depend on both the pH and redox potential compared to Heme I propionate D (plot b.1). Comparing the protonation probabilities of heme propionates A of Hemes II-V, Heme IV and Heme III propionates show strong dependence on both the pH and redox potential where as Heme II and V show dependence of the pH only. In heme propionates D of Hemes II-V, only Heme IV and V show strong dependence on pH and Heme II and III propionates are completely deprotonated and protonated respectively. The interaction between the heme charges depend on how deeply the hemes are buried in the protein as well as the distance between the hemes. The stabilization of the hemes by propionic acids can 4.4. Conclusions 99 Figure 4.9. Interaction energies of Hemes (I-V) with heme propionates A and D in cytochrome cnitrite reductase. The plot shows the interaction energies of Hemes (I-V) with heme propionates A and D in cytochrome cnitrite reductase. also have influence on the redox potential of the hemes. The interaction energies between the hemes are higher compared to the heme propionates. The interaction energies between Heme III and IV is around 3 kcal/mol due to the staking position of Heme III and IV with respect to Heme I and V. The analysis of the interaction between hemes show that they have influence on the heme mid-point potentials. The coupling of the pKashift of propionates with redox potentials of the hemes are supported from the structure because of the proximity of hemes to these propionates. 4.4 CONCLUSIONS Electrostatic calculations were performed on a simpler, better resolved crystal structure of cytochrome cnitrite reductase before doing calculations on the complex system cytochrome c oxidase. The redox potentials were calculated and in the range between –270 to –30 mV. The redox potentials are almost same for the oxidized and reduced enzyme except for Heme III. The oxidation and protonation probabilities of the hemes and heme propionates were calculated. The redox potentials of the hemes are sensitive to their location and are also influenced by the protonation probabilities of the surrounding titratable groups and as well as the propionates in the protein. Heme I, III and IV are close enough to allow direct π-electron interactions of the porphyrin ring. The mid-point potential of the Hemes II-V decreases with increase in pH. The oxidation probabilities of the Hemes II-V show less dependence on pH and Heme III and IV shows strong dependence. The oxidation probabilities show that the redox states of the hemes are strongly 100 Protonation and redox potentials of cytochrome cnitrite reductase depend on the pH. The protonation probabilities of Heme I propionates are strongly depend on the redox states of other hemes. On comparing the protonation probabilities of propionates A of Hemes II-V, propionates of Heme IV and Heme III show strong dependence on both the pH and redox potential and where as Heme II and V show dependence on the pH only. In heme propionates D of Hemes II-V, only Heme IV and V show strong dependence on pH and Heme II and III propionates are completely deprotonated and protonated respectively. The interaction between the hemes and heme propionates play a major role in tuning the redox potentials of the hemes. 4.4. Conclusions 101 a.2. a.3. a.4. a.5. b.2. b.3. b.4. b.5. Figure 4.10. Protonation probabilities of the heme propionates A and D of Hemes II-V. (a.2)-(a.5) Calculated protonation probabilities of heme propionates A of Hemes II-V. (b.2)-(b.5) Calculated protonation probabilities of heme propionates D of Hemes II-V. The probabilities are color coded as indicated by the scale next to the diagram. The green color indicates that the protonation probability of the residues is about 0.5. 102 Protonation and redox potentials of cytochrome cnitrite reductase CHAPTER 5 PKaCALCULATIONS OF BINUCLEAR CENTER OF CYTOCHROME cOXIDASE USING DFT CALCULATIONS Knowledge of the pKavalues of ionizable groups are important for an understanding of many phenomena of chemistry of both the gas phase and the solution. The pKavalues are of particular interest in elucidating the reaction mechanisms, especially those involving proton transfers and interpreting the binding of substrates or inhibitors to enzymes. Experimental determinations of individual pKavalues are difficult in complex systems. One case in point concerns the direct measurement of pKavalues of titrating groups of catalytically important residue or substrate in the enzyme-substrate complexes [101]. Cytochrome coxidase contains four redox centers CuA, heme a, heme a3and CuB. It catalyzes the reduction of oxygen to water. In the course of its catalytic activity, electrons delivered by cytochrome care transferred via CuAand heme ato the Fea3-CuBthe binuclear center, where the reduction of oxygen takes place. For every molecular oxygen a total of eight protons are consumed from the matrix side of the mitochondrion: four protons are delivered to the binuclear center for water formation and another four protons are translocated across the membrane (see Eq. (5.1)). 4 cyt cred + 8 H+ (in)+ O2−→ 4 cyt cox + 2 H2O + 4 H+ (out)(5.1) A detailed picture of the underlying mechanism of how the protons are translocated to the binuclear center and how they reach the opposite side of the enzyme is a matter of debate [57, 179–181]. The analysis of the many available crystallographic structures of cytochrome c oxidase does not show any specific pathways connecting the D-channel to the binuclear center or leading protons to the exterior side of the membrane. Recent studies have proposed the direct involvement of the one of the CuBhistidine ligands His291 in bovine heart or His334 in R. sphaeroides, a proton loading site in the pumping mechanism [1, 2, 8, 12]. Based on electrostatic calculations Stuchebrukhov et al. proposed that CuBbound imidazole ring of His334 play a key role in the proton pumping mechanism [1, 2]. The central feature of the proposed mechanism is that the pKavalues of the imidazole vary significantly depending on the redox state of the metals in the binuclear center [5, 6]. The proposed model states 103 104 pKacalculations of binuclear center of cytochrome coxidase using DFT calculations that upon reduction of the binuclear center, His334 deprotonates at its Nδ1[5]. From DFT calculations on the models of CuBin aqueous solution, these authors found that the pKavalue of His334 is 8.6 for the oxidized CuBcenter and 13.2 for the reduced CuBcenter. Another DFT study [182] on the CuBcenter obtained a pKavalue of above 13.5 for His334 in the oxidized state making it unlikely that a redox-coupled protonation state change of His334 directly involved in the proton pumping mechanism. The discrepancy between these two studies lead to some debate in the field [4, 183]. However in the two studies [5, 182], the crystal coordinates were not fully relaxed leading to unrealistic bond lengths. Moreover in the later study [183], it was not considered that some residues of the cytochrome coxidase can adopt a protonation that deviates from their usual protonation in aqueous solution. Siegbahn and co-workers, performed hybrid density functional theory calculations using B3LYP functional to study the energetics of the proton translocation in cytochrome coxidase [184]. They also computed the redox potentials of the metal centers and the tyrosyl radicals, as well as pKa values of important groups along the translocation paths [184, 185]. Ryde and co-workers systematically investigated how the axial ligand in heme proteins influences the geometry, electronic structure and spin states of the active site using density functional B3LYP method and medium-sized basis sets [186]. To investigate the reaction mechanism and as well as the role of the histidines bound with the CuBcenter and the histidine coordinating the heme a3, we calculated the pKavalues of the Cu-bound imidazole in various CuBand Fea3complexes using density functional theory in combination with classical continuum electrostatics models. In the previous studies only, the pKavalue of His334 was determined [5, 6, 182] and not all possible microscopic pKa values were calculated in order to check for internal consistency. In the present study, we calculated all microscopic pKavalues of the CuBand heme a3center both in the reduced and in the oxidized states [187]. The Finite Difference Poisson Boltzmann (FDPB) method and conductor-like polarizable continuum model (C-PCM) are used to determine the solvation free energies in aqueous solution. The dependence of the solvation free energies on the parameters are explored by comparing these two solvation models. The present study attempts to find reliable method to obtain proper pKavalues of the CuBcenter in aqueous solution. 5.1 PKaCALCULATIONS pKacalculations are used mainly in studies of enzyme mechanisms. Brief introduction of acid-base equilibria reactions are given in section 2.1.1 and 2.1.2 [62]. It is required to have reliable and accurate means of calculating relative and/or absolute pKavalues and proper understanding of the factors involved in pKacalculations. The pKavalues are very sensitive to the charges used for the solvation free energy calculations. In this chapter, the pKacalculations on CuBand heme a3center in the aqueous solution are reported. The pKa calculations in cytochrome coxidase are reported and discussed in chapter 6. The accuracy of the pKacalculations depend on different factors such as the solvation energies of the reactant and product, the solvation energy used for the proton and the model used for the pKa calculations. The calculation of absolute pKavalues of species in solution is a delicate task, which necessitates the use of high-level QM methods capable of attaining chemical accuracy [63]. The treatment of solvation represents the least reliable aspects of pKacalculations [182]. The heart of the pKaestimation methods are the calculated atomic charges in case of FDPB 5.2. Computational Methods 105 method. To gauge the dependence of pKavalues on the level of theory used to describe the solvation, FDPB method and C-PCM [121, 188] are used to determine the solvation free energies in bulk water. 5.2 COMPUTATIONAL METHODS 5.2.1 DENSITY FUNCTIONAL CALCULATIONS The starting structures of the CuBand heme a3centers were taken from the X-ray crystal structure of bovine heart cytochrome coxidase obtained by Yoshikava et al. at 2.3 ˚ A resolution (2OCC in the protein data bank) [16]. The minimal model required to calculate the pKaof a CuBcenter consists of the Cu ion, methylimidazole representing the coordinating histidines 284, 333 and 334, a methyl group representing Tyr288 (which is cross-linked to His284) and the bound H2O molecule as CuBligand. The methyl groups attached to the imidazole rings represent the Cβof histidines. A full description of Tyr288 is not included because the calculations are performed only to obtained the pKacalculations of histidines and H2O coordinating to CuBcenter. The model compound considered for the CuBcenter is given in Figure 5.1. The model for the heme a3center consists of the heme a3and an axial histidine which is modeled by methylimidazole. The sixth coordination is modeled with H2O molecule as in the CuBcenter. The hydrophobic hydroxyethyl-farnesyl group is truncated next to the hydroxyl group. The model compound considered for the heme a3center is given in Figure 5.2. CuB H284 Y288 +2/+1 H333 H334 Figure 5.1. The model compound of CuBcenter. The histidines coordinating the copper ion are modeled by methyl imidazoles. The cross-linked tyrosine 288 is replaced by methyl group. The fourth coordination is the H2O molecule. The proton binding sites are indicated by white spheres. The CuBcenter has three proton binding sites which are the H2O molecule, Nδ1of His333, and Nδ1of His334. The three proton binding sites lead to eight possible microscopic states and twelve proton equilibrium reactions [73]. In heme a3center, there are two proton binding sites the H2O molecule and the Nδ1of His419. The two proton binding sites lead to four possible microscopic states and four proton equilibrium reactions [73]. All possible protonation states 112 pKacalculations of binuclear center of cytochrome coxidase using DFT calculations Step Rea. State Prod. State TZP TZ2P+axc-BLYPbunconstraintc ∆∆Gsolv pKa∆∆Gsolv pKa∆∆Gsolv pKa∆∆Gsolv pKa 1Cu+2 B-111 Cu+2 B-011 98.77 6.9 99.00 7.9 97.56 6.5 80.93 6.5 Cu+2 B-111 Cu+2 B-101 96.19 16.1 96.76 17.6 96.38 13.4 80.31 15.8 Cu+2 B-111 Cu+2 B-110 97.60 15.1 96.02 14.9 96.33 11.5 80.11 13.3 2a Cu+2 B-110 Cu+2 B-010 17.47 6.4 22.21 7.3 19.32 7.6 24.10 8.4 Cu+2 B-101 Cu+2 B-001 21.45 4.4 19.22 5.5 21.55 5.0 25.13 3.3 2b Cu+2 B-110 Cu+2 B-100 20.42 10.9 23.64 13.9 21.09 10.0 27.71 12.9 Cu+2 B-011 Cu+2 B-001 18.87 13.6 16.98 15.2 20.37 11.9 24.52 12.7 2c Cu+2 B-011 Cu+2 B-010 16.30 14.6 19.23 14.4 18.09 12.6 23.28 15.2 Cu+2 B-101 Cu+2 B-100 21.83 9.9 22.90 11.2 21.04 8.1 27.50 10.3 3Cu+2 B-100 Cu+2 B-000 -48.81 4.6 -49.54 4.4 -47.95 5.6 -26.60 3.9 Cu+2 B-010 Cu+2 B-000 -45.86 9.1 -48.11 11.0 -46.18 8.0 -23.00 8.4 Cu+2 B-001 Cu+2 B-000 -48.43 10.1 -45.86 10.2 -48.46 8.7 -24.23 11.0 Table 5.2. The solvation energy difference between deprotonated and protonated forms and the pKavalues of all possible protonation equlibrium reactions in the oxidized state of CuBcenter. The solvation energy difference between deprotonated and protonated (∆∆Gsolv) forms and the pKavalues of all possible proton equlibrium reactions are given. All solvation energies were obtained by FDPB method using CHELPG charges. The solvation energies are given in kcal/mol. aPW91 with TZ2P+ basis set (diffusion functions for iron and copper). bPW91 with BLYP exchage-correlation functionals and TZP basis set. cPW91 with TZP basis set where the geometry is completely relaxed. His333 and His334. The second proton is removed from H2O in ‘step 2a’ and in ‘step 2b’ from His333 and from His334 in ‘step 2c’. In ‘Step 3’ the third proton is removed from H2O, His333 and His334. All the Tables reported in this section follow the same notation. The nomenclature used are shown in Figures 5.3, 5.4, 5.5 and 5.6 respectively. The pKavalues of all possible proton equilibrium reactions are calculated to check the internal consistency i.e., the deprotonation of the same proton binding site becomes more difficult when less protons are bound to the center. The solvation energies were calculated with different charges to obtain appropriate model pKafor the CuBand heme a3center which is required to calculate the average pKavalues in the protein. 5.3.1 PKaVALUES OF THE CUBCENTER The enthalpy difference between the deprotonated and the protonated form (∆Hdepro vac ) were calculated for all microscopic states of CuBcenter in the oxidized and in the reduced state using PW91 method with TZP, TZ2P+ basis set and also using BLYP as exchange and correlation functionals. Table 5.2 lists the solvation energy difference (∆∆Gsolv) between deprotonated and the protonated forms. The pKavalues are calculated based on Eq. (2.14). The pKavalues are very sensitive to the basis set used and also with different exchange correlation function- 5.3. Results and Discussion 113 Step Rea. State Prod. State PW91/TZP MDCmMDCdMDCq ∆∆Gsolv pKa∆∆Gsolv pKa∆∆Gsolv pKa 1Cu+2 B-111 Cu+2 B-011 112.57 17.1 97.55 6.1 97.40 5.9 Cu+2 B-111 Cu+2 B-101 117.00 31.4 95.40 15.5 93.70 14.3 Cu+2 B-111 Cu+2 B-110 115.00 27.9 95.60 13.7 96.60 14.4 2a Cu+2 B-110 Cu+2 B-010 23.85 11.1 20.49 8.6 15.54 5.0 Cu+2 B-101 Cu+2 B-001 29.28 10.2 22.14 5.0 20.30 3.6 2b Cu+2 B-110 Cu+2 B-100 36.45 22.6 22.79 12.6 17.40 8.7 Cu+2 B-011 Cu+2 B-001 33.71 24.5 19.99 14.4 16.60 11.9 2c Cu+2 B-011 Cu+2 B-010 26.28 21.9 18.54 16.2 14.74 13.4 Cu+2 B-101 Cu+2 B-100 34.45 19.1 22.99 10.7 20.30 8.8 3Cu+2 B-100 Cu+2 B-000 -46.06 6.6 -46.12 6.6 -49.58 4.1 Cu+2 B-010 Cu+2 B-000 -33.46 18.2 -43.82 10.6 -47.72 7.7 Cu+2 B-001 Cu+2 B-000 -40.89 15.6 -45.27 12.4 -49.58 9.2 Table 5.3. The solvation energy difference between deprotonated and protonated forms and the pKavalues of all possible proton equlibrium reactions in the oxidized state of CuBcenter. The solvation energy difference between protonated and deprotonated forms (∆∆Gsolv) and the pKavalues of all possible proton equlibrium reactions in the oxidized state of CuBcenter at PW91 method with TZP Basis set are given. The solvation energies are obtained from FDPB method using multipole derived charges. The solvation energies are given in kcal/mol. als. Comparing the pKavalues obtained from the different basis set and exchange functionals, the variation between the pKavalues are found to be irregular. The major problem observed in the pKavalues trend is that the pKavalues correspond to ‘step 3’ where a third proton is removed from the doubly deprotonated state is 4.6 for by TZP basis set and 4.4 for TZP with diffused functions (TZ2P+) and 5.6 with BLYP exchange and correlation functionals. The trends in the pKavalues show that the removal of proton from the double deprotonated states is favored which is not energetically plausible. The protonation equilibrium should show the correct trend i.e., the deprotonation of the same proton binding site becomes more difficult when less protons are bound to the center. On comparing the pKavalues of reactions which correspond to the removal of proton from H2O, Cu+2 B-111 →Cu+2 B-011, Cu+2 B-110 →Cu+2 B-010 and Cu+2 B-100 →Cu+2 B-000, the removal of third proton has the lower pKavalue compared to ‘step 1’ and ‘step 2a’. The solvation energy difference between the reactions are also differ very much. The solvation energy difference are very high for the reactions in ‘step 1’ where a single proton is removed from the fully protonated state. The total charges for the reactants are +2 and the product is +1 in ‘step 114 pKacalculations of binuclear center of cytochrome coxidase using DFT calculations Micro. State TZP TZ2Paxc-BLYPbUnconstraintcMDCmdMDCddMDCqd Cu+2 B-111 -150.80 -150.69 -150.17 -150.16 -175.60 -149.7 -154.90 Cu+2 B-011 -52.03 -51.69 -52.61 -53.90 -63.03 -52.15 -57.50 Cu+2 B-101 -54.61 -53.93 -53.79 -54.74 -58.60 -54.30 -61.20 Cu+2 B-110 -53.20 -54.67 -53.84 -55.02 -60.60 -54.10 -58.30 Cu+2 B-001 -33.16 -34.71 -32.24 -34.50 -29.32 -32.16 -40.90 Cu+2 B-010 -35.73 -32.46 -34.52 -36.18 -36.75 -33.61 -42.76 Cu+2 B-100 -32.78 -31.03 -32.75 -31.27 -24.15 -31.31 -40.90 Cu+2 B-000 -81.59 -80.57 -80.70 -81.51 -70.21 -77.43 -90.48 Table 5.4. Solvation energies of each microscopic state of CuBcomplex in oxidized state. The solvation energies of the microscopic state of CuBcomplex obtained by FDPB method. aPW91 with TZ2P+ basis set (diffusion functions for iron and copper). bPW91 with BLYP exchage-correlation functionals and TZP basis set. cPW91 with TZP basis set where the geometry is completely relaxed during optimization. dPW91 with TZP basis set and MDCs are used for solvation energy calculations. From a−cCHELPG charges are used for solvation enegy calculations. All the solvation energies are given in kcal/mol. 1’ reactions. Use of higher basis set and different exchange correlation functionals does not help very much to improved the trends in the pKavalues. To find the proper reason for such trends in the pKathe solvation energies were calculated with different charges. The charges obtained from the multipole moment are used instead of CHELPG charges. Table 5.3 show the solvation energy difference between deprotonated and protonated forms where the solvation energies were obtained from FDPB method using multipole derived charges (MDC) (see section 3.1.2). The pKavalues are very sensitive to the charges used. The solvation energies obtained from monopole derived charges are very high compared to the solvation energies obtained from dipole and quadrupole charges. Obviously, the pKavalues are very high with MDCm. On comparing the Tables 5.2 and 5.3, the pKavalues show almost the same trend. The solvation energies obtained from each microscopic states of the CuBcenter using FDPB methods with different basis set, exchange correlation functionals, different charges are given in Table 5.4. The solvation energies obtained from the the CHELPG charges are similar with ∼1-2 kcal/mol difference with different basis set and exchange correlation functionals. The solvation energies from the MDCmcharges are much higher compared to the solvation energies obtained from CHELPG charges. The pKavalues were calculated for the reduced state of the CuBcenter to obtain the pKavalues of H2O, His333 and His334. The gas phase free energies are obtained for all microscopic state of CuBcenter in reduced state using PW91/TZP basis set. The enthalpy difference between deprotonated and protonated form (∆Hdepro vac ), the solvation energies of the deprotonated (∆GA− solv) and protonated (∆GAH solv), and the calculated pKavalues are given in Table 5.5. The pKavalues of the reduced state of the CuBcenter also show the same trend like oxidized state of CuBcenter. The pKavalues corresponding to ‘step 2a’ are lower than the reactions in ‘step 1’. The solvation energies 5.3. Results and Discussion 115 Step Rea. State Prod. State PW91/TZP ∆GAH solv ∆GA− solv ∆Hdepro vac ∆∆Gsolv pKa 1Cu+1 B-111 Cu+1 B-011 -49.55 -34.35 -2.03 15.20 19.9 Cu+1 B-111 Cu+1 B-101 -49.55 -36.06 -8.59 13.59 14.7 Cu+1 B-111 Cu+1 B-110 -49.55 -33.54 -9.14 16.01 15.3 2a Cu+1 B-110 Cu+1 B-010 -34.35 -77.56 67.82 -49.60 9.0 Cu+1 B-101 Cu+1 B-001 -33.54 -78.32 47.80 -42.29 13.0 2b Cu+1 B-110 Cu+1 B-100 -77.56 -197.25 48.48 -44.78 27.5 Cu+1 B-011 Cu+1 B-001 -82.98 -197.25 41.24 -43.90 16.7 2c Cu+1 B-011 Cu+1 B-010 -78.32 -197.25 60.71 -48.79 28.1 Cu+1 B-101 Cu+1 B-100 -34.35 -82.98 47.93 -42.36 19.8 3Cu+1 B-100 Cu+1 B-000 -36.06 -77.56 143.20 -118.93 14.2 Cu+1 B-010 Cu+1 B-000 -33.54 -82.98 123.86 -114.11 24.4 Cu+1 B-001 Cu+1 B-000 -36.06 -78.32 143.33 -119.00 13.6 Table 5.5. The enthalphy (∆Hdepro vac ) and the solvation energy difference (∆∆Gsolv) between protonated and deprotonated forms, the solvation energies of the protonated (∆GAH solv) and deprotonated (∆GA− solv) form, and the calculated pKavalues of CuBcenter in the reduced state. The enthalphy of differendce (∆Hdepro vac ) between protonated and deprotonated form, the solvation energies of the protonated (∆GAH solv) and deprotonated (∆GA− solv) form, the solvation energy difference between deprotonated and protonated form (∆∆Gsolv), the calculated pKavalues of CuBcenter in the reduced stateform at DFT method with TVP basis set are given. The solvation energies are obtined by FDPB method using CHELPG charges. The enthalphies and solvation free energies are given in kcal/mol. are very low for reactions correspond to ‘step 3’. The ‘step 1’ reactions show positive solvation energies. The enthalpy difference between the deprotonated and the protonated forms are positive for the reactions corresponding to ‘step 3’ where the third proton is removed from the doubly deprotonated state. The results show that the pKavalues depend on the basis sets, different exchange and correlation functionals and also the mode of the optimization. But the lower pKavalues obtained for the triply deprotonated reactions may be because of the solvation model used for the solvation energy calculations. All the solvation energies from Table 5.2 to 5.4 were calculated using the FDPB method. Though this method was successful in many cases, [101–104] for the calculation of the pKa values of CuBcomplexes, this method seems to be insufficient. To have better understanding of the factors influencing the pKavalues, the pKacalculations are also performed using the B3LYP method with 6-31+G* and 6-31G* basis sets. The solvation energies were calculated by C-PCM and FDPB method using different charges like Mulliken, NPA, CHELPG, MK and RESP. The pKavalues obtained at B3LYP level of theory with 6-31+G* basis set and the solvation energy difference between deprotonated and protonated forms of the oxidized state of CuB 116 pKacalculations of binuclear center of cytochrome coxidase using DFT calculations Step Rea. State Prod. State B3LYP/6-31+G* CHELPG MK RESP ∆∆Gsolv pKa∆∆Gsolv pKa∆∆Gsolv pKa 1Cu+2 B-111 Cu+2 B-011 102.49 6.9 101.64 6.3 101.25 6.0 Cu+2 B-111 Cu+2 B-101 105.21 20.4 105.32 20.5 105.39 20.5 Cu+2 B-111 Cu+2 B-110 104.95 19.0 104.94 19.0 105.16 19.1 2a Cu+2 B-110 Cu+2 B-010 18.65 2.1 18.31 1.9 18.05 1.7 Cu+2 B-101 Cu+2 B-001 22.20 -0.5 20.14 -2.0 20.47 -1.8 2b Cu+2 B-011 Cu+2 B-001 24.92 13.0 23.82 12.2 24.61 12.7 2c Cu+2 B-011 Cu+2 B-010 21.11 14.2 21.61 14.5 21.96 14.8 3Cu+2 B-010 Cu+2 B-000 -37.33 9.5 -38.06 8.9 -38.24 8.8 Cu+2 B-001 Cu+2 B-000 -41.14 10.7 -40.27 11.3 -40.89 10.9 Table 5.6. The solvation energy difference between deprotonated and protonated forms (∆∆Gsolv) and pKavalues of all possible proton equlibrium reactions in the oxidized state of CuBcenter. The pKavalues are obtained at B3LYP/6-31+G* basis set. The solvation energies are calculated with CHELPG, MK and RESP charges. The solvation energies are calculated using FDPB method. The solvation energies are given in kcal/mol. center are given in Table 5.6. Initially the solvation energies were calculated with CHELPG, MK and RESP charges. The RESP charges are obtained by using hyperbolic restraint and ESP from MK method. The solvation energies were calculated using FDPB method. All the microscopic states of the CuBcenter were optimized with B3LYP/6-31+G* basis set. The ESP charges CHELPG and MK charges were obtained by fitting the ESP obtained from B3LYP/631+G* by CHELPG and MK method. The microscopic state Cu+2 B-100 was unstable and did not converge. The reactions involving Cu+2 B-100 microscopic state is not reported in the Table 5.6. On comparing the solvation energy difference (∆∆Gsolv) between the reaction and product states, the solvation energies correspond to ‘step 1’ are more positive compared to other steps. The pKavalues of ‘step 2a’ are very low and even negative for the reaction Cu+2 B-101 →Cu+2 B001. According to the pKavalues, the removal of second proton from the singly deprotonated state is favored compared to the removal of single proton from the fully protonated state. The pKavalues obtained using CHELPG and MK are almost similar and in some cases they vary up to 1 pKaunit. The pKavalues obtained from RESP charges are close to the pKavalues obtained from MK charges. Either the level of theory or the basis set with diffusion and polarization function does not help to improve the pKatrends. To check the effect of diffusion function on the pKacalculations, the calculations were done at B3LYP level of theory with 6-31G* basis set without diffusion functions. The solvation energy difference between deprotonated and protonated forms, the pKavalues of the oxidized state of CuBcenter are given in Table 5.7. The solvation energies were calculated with Mulliken, NPA, CHELPG and MK 5.3. Results and Discussion 117 Step Rea. State Prod. State B3LYP/6-31G* Mulliken NPA CHELPG MK ∆∆Gsolv pKa∆∆Gsolv pKa∆∆Gsolv pKa∆∆Gsolv pKa 1Cu+2 B-111 Cu+2 B-011 104.54 11.3 103.52 10.6 105.71 12.1 105.91 12.3 Cu+2 B-111 Cu+2 B-101 93.22 20.0 86.91 15.5 94.38 20.9 94.36 20.9 Cu+2 B-111 Cu+2 B-110 94.94 19.4 87.92 14.3 94.92 19.4 95.71 20.0 2a Cu+2 B-110 Cu+2 B-010 33.24 13.3 32.15 12.5 31.66 12.1 30.73 11.5 Cu+2 B-101 Cu+2 B-001 32.22 8.7 28.38 5.9 32.89 9.2 32.52 8.9 2b Cu+2 B-110 Cu+2 B-100 27.62 15.8 22.07 11.7 30.33 17.7 27.73 15.8 Cu+2 B-011 Cu+2 B-001 20.90 17.4 11.77 10.8 21.56 17.9 20.97 17.5 2c Cu+2 B-011 Cu+2 B-010 23.64 21.4 16.55 16.2 20.87 19.4 20.53 19.1 Cu+2 B-101 Cu+2 B-100 29.34 15.1 23.08 10.6 30.87 16.2 29.08 14.9 3Cu+2 B-100 Cu+2 B-000 -32.74 13.0 -35.64 10.9 -35.53 11.0 -33.71 12.3 Cu+2 B-010 Cu+2 B-000 -38.36 15.5 -45.72 10.2 -36.86 16.6 -36.71 16.7 Cu+2 B-001 Cu+2 B-000 -35.62 19.5 -40.94 15.7 -37.55 18.1 -37.15 18.4 Table 5.7. The solvation energy difference between deprotonated and protonated forms (∆∆Gsolv) and pKavalues of all possible proton equlibrium reactions in the oxidized state of CuBcenter The pKavalues are obtained at B3LYP/6-31G* basis set. The solvation energies are calculated with Mulliken, NPA, CHELPG and MK charges. The solvation energies are calculated using FDPB method. The solvation energies are given in kcal/mol. charges. The solvation energies were also calculated with the RESP charges with two different restraints hyperbolic and harmonic restraint using Mulliken charges as initial charges. The pKavalues obtained are much better compared to the pKavalues obtained at B3LYP/6-31+G* basis set. But the pKavalues do not show correct trend i.e., the deprotonation of the same proton binding site becomes more difficult when less protons are bound to the center. The pKavalues obtained using the CHELPG and MK charges are comparable. The solvation energy is more positive for the single deprotonated step compared to the others steps. The solvation energies are negative for the reactions where the third protons are removed from the doubly protonated state. The pKavalues obtained with either Mulliken or NPA charges show the same trend in the pKavalues as in pKavalues obtained from ESP based charges. The solvation energy difference between deprotonated and deprotonated forms, pKavalues of all possible proton equilibrium reactions in the oxidized state, and the solvation energies calculated with RESP charges are given in Table 5.8. The solvation energies were calculated using FDPB method. The RESP charges are obtained by using the Mulliken and NPA charges as initial charges for harmonic and hyperbolic restraint. The electrostatic potential from MK charges were used. The pKavalues obtained by using RESP charges also show the same trend as other ESP charge methods. To get the proper understanding of the influence of the solvation energy on the pKacalculations, the C-PCM was used to calculate the solvation free 118 pKacalculations of binuclear center of cytochrome coxidase using DFT calculations Step Rea. State Prod. State B3LYP/6-31G* RESPaRESPb qwt=0.005 qwt=0.01 qwt=0.005 qwt=0.01 ∆∆Gsolv pKa∆∆Gsolv pKa∆∆Gsolv pKa∆∆Gsolv pKa 1Cu+2 B-111 Cu+2 B-011 104.91 11.6 104.59 11.3 105.83 12.2 105.64 12.1 Cu+2 B-111 Cu+2 B-101 94.85 21.2 94.23 20.7 94.50 21.0 94.52 21.0 Cu+2 B-111 Cu+2 B-110 96.89 20.8 96.76 20.7 95.68 19.9 95.70 20.0 2a Cu+2 B-110 Cu+2 B-010 30.06 11.0 30.36 11.2 31.02 11.7 31.12 11.7 Cu+2 B-101 Cu+2 B-001 33.14 9.3 33.38 9.5 33.37 9.5 33.81 9.8 2b Cu+2 B-110 Cu+2 B-100 29.36 17.0 29.65 17.2 28.95 16.7 29.16 16.9 Cu+2 B-011 Cu+2 B-001 23.08 19.0 23.02 18.9 22.04 18.2 22.69 18.7 2c Cu+2 B-011 Cu+2 B-010 22.04 20.2 22.53 20.6 20.87 19.4 21.18 19.6 Cu+2 B-101 Cu+2 B-100 31.40 16.6 32.18 17.2 30.13 15.7 30.34 15.8 3Cu+2 B-100 Cu+2 B-000 -36.33 10.4 -36.60 10.3 -34.31 11.9 -34.46 11.8 Cu+2 B-010 Cu+2 B-000 -37.03 16.5 -37.31 16.3 -36.38 17.0 -36.42 16.9 Cu+2 B-001 Cu+2 B-000 -38.07 17.7 -37.80 17.9 -37.55 18.1 -37.93 17.8 Table 5.8. The solvation energy difference between deprotonated and protonated forms (∆∆Gsolv) and pKavalues of all possible proton equlibrium reactions in the oxidized state of CuBcenter. The pKavalues are obtained at B3LYP/6-31G* basis set. The solvation energies are calculated with RESP charges. The solvation energies are calculated using FDPB method. The solvation energies are given in kcal/mol ahyperbolic restraint witih initial mulliken charges and electrostatic potential from MK method. bharmonic restraint with mulliken charges as initial charges and electrostatic potential from MK method. Restraint factor (qwt) of 0.005 and 0.01 are used for RESP calculations. energies. The solvation energy difference between deprotonated and protonated forms, pKa values obtained at B3LYP with TZVP, 6-31+G* and 6-31G* basis set are given in Table 5.9. The pKacalculations with C-PCM at B3LYP level with TZVP and 6-31+G* also show lower pKa values for the deprotonation reaction which are not energetically plausible. The pKavalues obtained at B3LYP with 6-31G* show the proper trend. The protonation equilibrium show the correct trend, i.e., the deprotonation of the same proton binding site becomes more difficult when less protons are bound to the center. The enthalpy of difference (∆Hdepro vac ), the solvation energy difference (∆∆Gsolv) and the calculated pKavalues of CuBcenter in the reduced state at B3LYP level with 6-31G* are given in Table 5.10. The solvation energies are obtained by using C-PCM. The pKain the reduced CuBcenter also show correct trend. The solvation energies obtained for all the eight microscopic states are given in Table 5.12 for comparison. The pKavalues obtained from B3LYP/6-31G* level of theory show that the solvation energies should be calculated in a self-consistency manner as in C-PCM. The FDPB model treats the interactions in a simpler way. The solute charge distribution is considered as classical entity and the electronic polarization due to solvent effects is discarded. The FDPB model suffers 5.3. Results and Discussion 119 Step Rea. State Prod. State C-PCM B3LYP/TZVP B3LYP/6-31+G* B3LYP/6-31G* ∆∆Gsolv pKa∆∆Gsolv pKa∆∆Gsolv pKa 1Cu+2 B-111 Cu+2 B-011 98.47 7.0 97.20 3.1 98.55 7.4 Cu+2 B-111 Cu+2 B-101 96.45 16.4 96.85 14.3 86.72 15.9 Cu+2 B-111 Cu+2 B-110 95.96 14.9 96.24 12.7 90.25 15.2 2a Cu+2 B-110 Cu+2 B-010 23.96 8.9 27.01 8.1 29.96 13.0 Cu+2 B-101 Cu+2 B-001 27.60 6.0 26.52 2.6 35.84 10.6 2b Cu+2 B-110 Cu+2 B-100 -a- -a- -a- -a27.84 18.1 Cu+2 B-011 Cu+2 B-001 25.58 15.4 26.17 13.9 24.01 19.1 2c Cu+2 B-011 Cu+2 B-010 21.45 16.8 26.05 17.7 21.66 20.8 Cu+2 B-101 Cu+2 B-100 -a- -a- -a- -a31.37 17.4 3Cu+2 B-100 Cu+2 B-000 -a- -a- -a- -a- -33.22 14.3 Cu+2 B-010 Cu+2 B-000 -32.48 16.3 -35.40 10.9 -35.34 19.4 Cu+2 B-001 Cu+2 B-000 -36.61 17.7 -35.52 14.8 -37.69 21.0 Table 5.9. The solvation energy difference between deprotonated and protonated forms (∆∆Gsolv) and pKavalues of all possible proton equlibrium reactions in the oxidized state of CuBcenter. The pKavalues are obtained at B3LYP with TZVP, 6-31+G* and 6-31G* basis set. The solvation energies are calculated with C-PCM. aThe structure of Cu+2 B-100 was unstable and did not converge at the given level of theory and the basis set used. by severe limitations to obtain the solvation energy of the CuBcomplexes. The present results show that in the solvation energy calculations the electronic polarization due to the solvent should be included explicitly. In C-PCM, the interaction energies between the solvent and the solute are included in the Hamiltonian itself. In C-PCM the electron distribution of the solute is obtained with the fixed nuclei in the presence of the reaction field of the solvent. Using the C-PCM for the oxidized state of the CuBmodel in the aqueous solution, a pKavalue of 15.9 and 15.2 are obtained for the first deprotonation reaction of His333 and His334 and the pKavalue of 7.4 for the H2O coordinating the copper. Obviously, the pKavalues of the two histidines are too high to allow their deprotonation in aqueous solution. The pKavalues of the reduced state is even higher because the negative charge of the electron is stabilizing the proton. The pKatrend is obtained correctly only in the combination of B3LYP/6-31G*/C-PCM level. 120 pKacalculations of binuclear center of cytochrome coxidase using DFT calculations Step Rea. State Prod. State B3LYP/6-31G* ∆GAH solv ∆GA− solv ∆Hdepro vac ∆∆Gsolv pKa 1Cu+1 B-111 Cu+1 B-011 -36.17 -15.89 304.45 20.28 32.5 Cu+1 B-111 Cu+1 B-101 -36.17 -12.90 280.73 23.27 17.7 Cu+1 B-111 Cu+1 B-110 -36.17 -12.37 277.77 23.80 15.9 2a Cu+1 B-110 Cu+1 B-010 -a- -a- -a- -a- -aCu+1 B-101 Cu+1 B-001 -12.90 -63.16 379.80 -50.26 36.2 2b Cu+1 B-110 Cu+1 B-100 -12.37 -47.50 344.09 -35.13 21.3 Cu+1 B-011 Cu+1 B-001 -15.89 -63.16 356.08 -47.27 21.1 2c Cu+1 B-011 Cu+1 B-010 -a- -a- -a- -a- -aCu+1 B-101 Cu+1 B-100 -12.90 -47.50 341.13 -34.60 19.5 3Cu+1 B-100 Cu+1 B-000 -47.50 -162.00 460.90 -114.50 48.4 Cu+1 B-010 Cu+1 B-000 -a- -a- -a- -a- -aCu+1 B-001 Cu+1 B-000 -63.16 -162.00 422.23 -98.84 31.7 Table 5.10. The enthalphy (∆Hdepro vac ) and the solvation energy difference (∆∆Gsolv) between protonated and deprotonated form and the calculated pKavalues of CuBcenter in the reduced state. The enthalphy difference ∆Hdepro vac between protonated and deprotonated, solvation energy difference between deprotonated and protonated form ∆∆Gsolv and the calculated pKavalues of CuBcenter in the reduced state at B3LYP/6-31G* level are given. The solvation energies are obtined by using C-PCM. aThe structure of Cu+1 B-010 was unstable and did not converge at the given level of theory and the basis set used. 5.3.2 PKaVALUES OF THE HEME a3CENTER The enthalpy of difference ∆Hdepro vac between deprotonated (HAH vac ) and protonated (HA− vac ), the solvation energies of the deprotonated (∆GA− solv) and protonated (∆GAH solv), solvation energy difference between deprotonated and protonated form ∆∆Gsolv and the calculated pKavalues of CuBcenter in the reduced and oxidized form of heme a3center are given in Table 5.11. The solvation energies were calculated using CHELPG charges in case of PW91 method and C-PCM solvation model was used in B3LYP/6-31G* level. The pKavalues obtained at B3LYP/6-31G* level are higher than the pKavalues obtained at PW91/TZP basis set. The pKavalues 19.3 and 10.7 correspond to the deprotonation of His419 and H2O respectively. The pKavalue of the doubly deprotonated and reduced heme a3center are even too high. The pKavalues of the histidine are too high to allow their deprotonation in aqueous solution. 5.4. Conclusions 121 Rea. State Prod. State PW91/TZPaB3LYP/6-31G* ∆Hdepro vac ∆∆Gsolv pKa∆Hdepro vac ∆∆Gsolv pKa Fea3+3-11 Fea3+3-01 268.65 17.33 7.8 268.65 21.47 10.7 Fea3+3-11 Fea3+3-10 277.26 19.7 10.4 277.26 24.56 19.3 Fea3+3-01 Fea3+3-00 343.72 -43 16.5 343.72 -33.52 25.4 Fea3+3-10 Fea3+3-00 335.11 -45.37 13.8 335.11 -36.61 16.8 Fea3+2-11 Fea3+2-01 364.65 -43.78 18.3 364.65 -45.85 31.7 Fea3+2-11 Fea3+2-10 349.27 -38.39 13.5 349.27 -36.00 27.7 Fea3+2-01 Fea3+2-00 402.75 -93.08 19.9 402.75 -79.27 35.2 Fea3+2-10 Fea3+2-00 418.13 -98.47 24.7 418.13 -89.12 39.3 Table 5.11. The enthalphy (∆Hdepro vac ) and the solvation energy difference between deprotonated and protonated form ∆∆Gsolv and the calculated pKavalues of CuBcenter in the reduced and oxidized states of heme a3center. The pKavalues are calculated using PW91/TZP and B3LYP/6-31G* basis sets. The solvation energies are calculated using CHELPG charges in case of PW91 method and C-PCM solvation model was used in B3LYP/6-31G* level. The solvation energies and enthalpies are given in kcal/mol. 5.4 CONCLUSIONS In the present study, we calculated all microscopic pKavalues of the CuBcenter combining PW91 and B3LYP with both FDPB and C-PCM solvation models, with the purpose of assessing the validity of the recently proposed pumping mechanism of cytochrome coxidase, based on the deprotonation of the His334 [1, 2]. The pKavalues of H2O, His333 and His334 were calculated in the aqueous solution to find the possibility of these ligands role in the proton pumping mechanism of cytochrome coxidase. Extensive studies were done to understand the influence of various factors on the pKavalues of the CuBcenter. The pKacalculations performed by two different density functional methods PW91 and B3LYP and solvation models (FDPB and C-PCM) show that in solvation energy calculations the electronic polarization due to the solvent should be included explicitly, which is lacking in FDPB method. Due to this limitations, the pKavalues calculated using the solvation energies obtained from FDPB method do not show the trend correctly i.e., the removal of the third proton is favored compared to other reactions [193, 194] which is energetically not plausible. The B3LYP/6-31G*/C-PCM level should be used to get correct pKatrend for CuBwhere solvation energies were calculated by C-PCM in which the electronic polarization of the solute due to the solvent reaction field is obtain by self-consistent manner. The experimental pKavalue for the deprotonation of imidazole is approximately 14 [193, 194]. Our calculations, in aqueous solution on small CuBmodels do not show any significant change of the pKadue to the ligation of the imidazoles to Cu ion. Further more, the results obtained from reduced CuBcenter revels a small increase in the pKavalues of imidazoles in aqueous solution. The pKavalues of 15.9 and 15.2 for the first deprotonation reaction of His333 and 128 pKacalculations of CuBligands in cytochrome coxidase In the transition between the F→Fpstates, a chemical proton (proton involved in oxygen reduction) enters the binuclear center and converts the hydroxyl CuBligand into water thereby increasing the total charge of the site from +1 to +2. This state of the catalytic cycle is supposedly connected to the first deprotonation of the His334 imidazole ring. The coupled electron-proton transfer reaction leads to H0state. It was proposed [1] that H0is likely in dynamic equilibrium with the H state in which iron is in the +3 state and copper ion is in the +2 state and each of the metals is ligated to the hydroxyl group. The transition between the H and O states is connected to a second deprotonation of His334. The protonation of the hydroxyl group bound to heme a3leads to the formation of O state, which is associated with the pumping event. The O state can evolve in two different ways. The proton entering the binuclear center can protonate the hydroxyl group of CuBcenter which is energetically least favorable (O∼state). The another possible pathway (OH) results in the formation of E state. The E state was proposed [1] to be thermodynamically unstable and converted to Ep. In E state, the redox state of heme a3changes from +3 to +2. In Epstate, the chemical proton enters the active site and converts the hydroxyl group to water, triggering the third deprotonation of His334. The next incoming proton reloads the His334 site, closing the catalytic cycle with removal of water molecule. In the two studies [1–3, 5, 6, 182], the crystal coordinates of the CuBcenter was not fully relaxed during the geometry optimizations which can lead to unrealistic bond lengths. Moreover in later studies [182], it was not considered that some residues of the cytochrome coxidase can adopt a protonation deviating from their usual protonation in aqueous solution. In the two studies, only the deprotonation of His334 was determined and not all possible microscopic proton equilibrium reactions in the protein are considered. There is a possibility that His333, another ligand of CuBcenter can serve as a proton loading site for subsequent pumping which was not examined. To understand the coupling of the electron transfer and proton translocation, DFT and continuum electrostatic calculations were performed on R. sphaeroides cytochrome coxidase. Initially all the microscopic pKavalues of the bound H2O, His333 and His334 in CuBcenter were calculated in aqueous solution by combining DFT calculations with PCM model (see chapter 5). In order to take the influence of the protein environment into account, the average pKavalues of H2O, His333 and His334 ligands bound to CuBwere calculated in the protein by considering the protonation of the protein at pH=7. The CuBcenter was considered in both reduced and oxidized state. The heme a3center was considered with all possible ligand states and ligand unbound states. The electrostatic calculations were performed by solving the Poisson Boltzmann Equation by finite difference method. 6.2 STRUCTURE PREPARATION AND MODELS Over the past years, substantial progress has been made in the understanding of the structure and function of cytochrome coxidase. A landmark in the field of cytochrome coxidase research was the determination of the three-dimensional structures of the bacterial cytochrome coxidase from the soil bacterium P. denitrificans [12, 14] and its mammalian counterpart from 6.2. Structure preparation and models 129 bovine heart mitochondria [16–18]. The first crystals of the bacterial enzyme were obtained using the four subunit cytochrome coxidase [12]. Later, an improved crystal form of the functionally active two-subunit cytochrome coxidase complex with the Fvantibody fragment as the four-subunit cytochrome coxidase was observed and the structure was determined at 2.7 ˚ A resolution [14]. The structures of the bacterial and the mitochondrial enzymes are surprisingly similar. The core parts (subunits I, II, and III) of the two crystal structures look nearly identical at the atomic level. The following section deals with the structures and necessary preparations for the electrostatic calculations of cytochrome coxidase. 6.2.1 PREPARATION OF X-RAY STRUCTURE OF PROTEIN The electrostatic calculations were performed on the X-ray structure of R. sphaeroides (PDB ID:1M56). The structure of cytochrome coxidase have been determined at 2.8 ˚ A resolution by Iwata et al. [15]. The crystal structures were determined for wild type and mutant by the replacement of glutamate-286 of subunit I by glutamine. The electrostatic calculations were performed on the wild type enzyme. The overall structure of the subunit I-IV of R. sphaeroides is very similar to those of the P. denitrificans while subunit I-III are similar to the corresponding subunits of the bovine heart cytochrome coxidase. The environments around the redox-metal centers (the binuclear center, heme aand CuA) and the non-redoxactive centers (Mg2+ and Ca2+) are structurally very similar to those of both the bovine heart and P. denitrificans. Six phospholipid molecules were identified in the X-ray structure of R. sphaeroides and are assigned as phosphatidylethanolamines from the shapes of the electron density. These phospholipid molecules are also included in the electrostatic calculations. The electrostatic calculations were performed in the monomer of cytochrome coxidase. Structures were prepared for electrostatic calculations using the CHARMM [153] molecular modeling package. Hydrogen atom positions were generated using the HBUILD algorithm implemented in the CHARMM program. All the heavy atoms were fixed and energy was minimized using the CHARMM forcefield. The energy minimizations were performed by 500 steepest decent (SD) steps, followed by 500 conjugate gradient (CG) steps. Partial atomic charges for atoms of the standard amino acids were taken from the CHARMM parameter set. The partial atomic charges for the redox-metal centers (the binuclear center, heme aand CuA) and the non-redox-active centers (Mg2+ and Ca2+) were obtained from quantum chemical calculations (see section 6.3). The protein structures were minimized with the crystallographic water molecules. In electrostatic calculations the crystallographic water molecules were removed since the construction of water hydrogens would arbitrarily assign a certain orientation to the water molecules that would affect electrostatic calculations. The C-terminus of subunits I-IV of cytochrome coxidase were not resolved and a neutral blocking group (acetly group) was added to the C-terminus. CONSTRUCTION OF A MEMBRANE MODEL. In order to include the effect of the membrane environment into the electrostatic calculations, the protein complex was embedded into a cylindrical shaped belt of uncharged dummy atoms where the dummy atoms model the hydrophobic membrane core (see Figure 6.2). The dummy atoms which mimics the membrane environ- 130 pKacalculations of CuBligands in cytochrome coxidase ment were assigned a low dielectric constant in the electrostatic calculations. The Bondi radii were used for the protein atoms [175]. The cavities within the protein that show no connection to the membrane environment were not be filled with the dummy atoms. These cavities are considered to contain distorted water molecules and assigned a high dielectric constant. The structure of cytochrome coxidase with membrane was obtained from the orientations of protein in membranes (OPM) database [205–207]. The OPM database currently includes all unique structures of transmembrane protein complexes and selected monotopic, peripheral proteins and membrane-bound peptides from PDB with their calculated membrane boundaries. Coordinate files of the proteins with calculated membrane boundaries are available for download. The coordinates of cytochrome coxidase with calculated membrane boundaries were superimposed with the crystal structure of R. sphaeroides using the program superimpose. The membrane was constructed around this superimposed structures. The Vol program [208] was used to construct a membrane around the protein. The thickness of the membrane model is 30 ˚ A . The radius of the contracted membrane is 50 ˚ A . A probe radius of 1.4 ˚ A was used for generating the protein surface. The membrane model was constructed for all structures considered in the electrostatic calculations. (A) (B) Figure 6.2. Membrane model. (A): View of cytochrome coxidase from R. sphaeroides along the membrane plane. The dummy atoms which mimics the membrane environment were assigned a low dielectric constant. (B): View from the cytoplasm. 6.2.2 REDOX CENTER MODELS In the reaction mechanism of cytochrome coxidase, four redox centers are involved. The reduction of molecular oxygen to water takes place in the binuclear center where the electrons and protons are delivered for oxygen reduction and four protons are translocated across the inner-mitochondrial membrane, a process which results in a membrane electrochemical pro- 6.2. Structure preparation and models 131 ton gradient. The proton translocation is coupled with the electron transfer which makes the reaction difficult to study. To incorporate the redox centers to the protein, the redox centers were considered for the quantum chemical calculations to get the charges with different oxidized and reduced states as well these models were used for the DFT calculations to obtain the gas phase and solvation energies. The models used for the redox centers are described below. DFT calculations were also performed in non-redox-active centers (Mg2+ and Ca2+) to obtain partial atomic charges. MODEL FOR CuACENTER. The CuAcenter contains two copper ions (see Figure 6.3) and are bridged by two cysteins sulfur atoms. The two copper ions of the CuAcenter are coordinated by two His, one Met, a backbone carbonyl oxygen of Glu, and two bridging Cys residues. The ligated cysteins were simplified as methyl thiolates. The histidines were modeled by the methyl imidazoles. The backbone carbonyl oxygen of Glu was included in the model compound. The amino acid side chains were cut at the Cαatom and their Cβatoms were fixed in their crystal structure positions. The initial coordinates were obtained from the crystal structure of R. sphaeroides. All the hydrogens were added using babel program. The CuAcenter model was used to obtain the charges for the calculations. The charges were obtained for both reduced and oxidized states of CuAcenter. CuA M263 C256 C252 H260 E254 H217 Figure 6.3. The model compound of CuAcenter The cysteins are simplified as methyl thiolates. The histidines are modeled by methyl imidazole. MODEL FOR HEME aHeme acenter is with two histidine residues as axial iron-ligands (see Figure 6.4). The histidine residues were modeled by methyl imidazole. The heme propionates were cut off and substituted by hydrogen atoms. The hydrophobic hydroxyethyl-farnesyl group was truncated next to the hydroxyl group. The calculations were performed both in the reduced and in the oxidized state. 132 pKacalculations of CuBligands in cytochrome coxidase heme a H421 H102 Figure 6.4. The model compound of heme a.The histidines coordinating the iron are modeled by methyl imidazole. The hydrophobic hydroxyethyl-farnesyl group was truncated next to the hydroxy group. The heme propionates were cut off and substituted by hydrogen atoms. MODEL FOR CuBCENTER. The model for the CuBcenter consists of Cu ion, methylimidazole model for coordinating histidines 284, 333, and 334 and the tyrosine 288 was modeled by methyl group (see Figure 6.5). The fourth coordination was modeled with H2O molecule and a hydroxyl group. The calculations were performed both in the reduced and the in the oxidized state. CuB H284 Y288 +2/+1 H333 H334 Figure 6.5. The model compound of CuBcenter. The histidines coordinating the copper ion are modeled by methyl imidazole. The cross linked tyrosine 288 is replaced by methyl group. The fourth coordination is the H2O molecule. The proton binding sites are indicated by white spheres. 6.3. Density Functional calculations 133 MODEL FOR HEME a3CENTER. The model for the heme a3center consists of the heme a3and a axial histidine which was modeled by methylimidazole (see Figure 6.6). The sixth coordination was modeled with H2O molecule (FeIII-H2O and FeII-H2O), hydroxyl group (FeIII-OH and FeII-OH), oxo-ferryl (FeIV=O) and ligand unbound state (FeIII and FeII). The hydrophobic hydroxyethyl-farnesyl group was truncated next to the hydroxyl group. The heme a3center was optimized both in oxidized and reduced states with all ligands mentioned above. His419 heme a3 Figure 6.6. The model compound of heme a3center. The histidines coordinating the iron are modeled by methyl imidazole. The hydrophobic hydroxyethyl-farnesyl group was truncated next to the hydroxyl group. The heme propionates were cut off and substituted by hydrogen atoms. The H2O molecule is considered in the six coordination position is shown. The proton binding sites are indicated by white spheres. 6.3 DENSITY FUNCTIONAL CALCULATIONS The DFT calculations [81, 172] were performed to obtain partial atomic charges of the redox centers CuA, heme a, heme a3and CuBrespectively. For redox centers CuA, heme a and heme a3the Perdrew-Wang 91 (PW91) calculations were performed with ADF 2004.01 [173] program. The partial atomic charges of CuBcenter were obtained by B3LYP method using 6-31G* basis sets using GAUSSIAN 03 program. The DFT calculations were performed on the oxidized states of CuAand heme a. Following seven states were considered for heme a3: aqua ferric state (FeIII-H2O, charge=+1, S=5/2), aqua ferrous state (FeII-H2O, charge=0, S=2), hydroxyl state (FeIII-OH, charge=0, S=5/2 and FeII-OH, charge=-1, S=2), oxo-ferryl state (FeIV=O, charge=0, S=2) and ligand unbound states (FeIII, charge=+1, S=5/2 and FeII, charge=0, S=2). The local density approximation (LDA) for exchange and correlation are based on the parametrization of Vosko, Wilk and Nausair [84]. The Perdrew-Wang 91 (PW91) [86] exchange and correlation functionals were used for the generalized gradient approximation (GGA). The numerical integration scheme used in this calculation was Voronoi polyhedron method developed by te Velde et al. with the accuracy parameter ACCINT set to its default value. A set of triple ζSlater type orbital (STO) was employed with single polarization function. The inner core shells were treated by the frozen core approximation. All the calculations were 134 pKacalculations of CuBligands in cytochrome coxidase done with a spin-unrestricted scheme. The optimization were performed by the quasi Newton method and the Hessian was updated with the Broyden-Fletcher-Goldfarb-Shanno strategy. The distributions of partial atomic point charges were computed by fitting the molecular electrostatic potentials calculated by ADF 2004.01 program [173]. The program chargefit was used to obtain point charges which is based on CHELPG algorithm [174]. The net charge of the molecule and the three Cartesian dipole moment components from PW91 calculations were adopted as constraints for the chargefit. Point charges were then computed by determining the electrostatic potential, how well the point charges effectively reproduce the electrostatic potential. The ESP charges were calculated on the cubic grid with uniform spacing of 0.2 and 3˚ A outer boundary around each atom of the molecule. The atoms were assigned the Bondi values of 1.7 for carbon, 1.2 for hydrogen, 1.55 for nitrogen, 1.5 for oxygen, 1.3 for iron and 1.8 for sulfur. To minimize the uncertainties in the fitting procedure the single value decomposition (SVD) method [209] was used to obtain a model with stable atomic charges. Partial atomic charges for the amino acids are taken from the CHARMM22 parameter set. The partial atomic charges for all the redox centers were derived from PW91 calculations. In addition to atomic coordinates and atomic charges, electrostatic calculations require radii. The radii used were taken from Bondi [175]. All microscopic states of CuBcenter were optimized with hybrid density functional calculations (B3LYP) using 6-31G* basis sets. The geometry optimizations were performed using GAUSSIAN 03 program. The single point calculations were performed on the optimized structures of the CuBcenter using B3LYP/6-31G*/C-PCM level. The electrostatic potential obtained from CPCM self-consistent reaction field were fitted by Merz-Kollman method to obtain the point charges. 6.4 ELECTROSTATIC CALCULATIONS The electrostatic calculations were performed on Fea3-CuBcomplexes (see Figure 6.7). The CuBcenter was allowed to adopt each of the eight states shown in Figure 5.3 and 5.4 respectively. Following seven possible states of heme a3center were considered: FeIV=O, FeIII-H2O, FeII-H2O, FeIII-OH, FeII-OH, FeIII and FeII. These seven states with different CuBcenter are shown in Figure 6.7. The CuAcenter and heme awere considered oxidized throughout the studies. Continuum electrostatic calculations were performed with the QMPB which solve the linear Poisson-Boltzmann equation by numerical finite difference method. The whole system is divided into three regions with three different dielectric constant of i= 1 for the active site (quantum region) in this case the CuBcenter, s= 80 for the solvent region, and p= 4 for the protein. Compared to the purely electrostatic dielectric constant = 2,p= 4 was adopted for the protein, which allows some mobility of the protein dipoles and accounts for some reorientational relaxation of the protein in an approximate way. The ionic strength was set to 0.1 mol/l and the absolute temperature is set to 300 K. The boundary between the interior and the exterior is defined as the solvent contact and re-entrant surfaces of 1.4 ˚ A spherical solvent probe rolling over the van der Waals surface of the protein. The atomic radii for the protein’s atom and the charges for the protein’s non-titrating atoms were taken from the polar hydrogen parameter set of CHARMM22. Electrostatic potential was calculated by first focusing the grid on the protein model compound and the second grid was centered on the titratable 6.4. Electrostatic calculations 135 Fe =O IV HH O HH O Fe −H O II Fe −H O III O H Fe −OH III O H Fe −OH II Fe III Fe II heme a3heme a3 heme a3 heme a3 CuBCuBCuB CuB CuB heme a3 heme a3 heme a3CuBCuB Cu His280 His333 His334 O HH 1+/2+ Fe Cu His280 His333 His334 O HH 1+/2+ Fe Cu His280 His333 His334 O HH 1+/2+ Fe Cu His280 His333 His334 O HH 1+/2+ Fe Cu His280 His333 His334 O HH 1+/2+ Fe Cu His280 His333 His334 O HH 1+/2+ Fe Cu His280 His333 His334 O HH 1+/2+ OFe 4+ 3+ 2+ 3+ 2+ 3+ 2+ 2 2 Figure 6.7. The Fea3-CuBcomplexes considered in the present study. The Fea3-CuB complexes considered for the present study are shown. The nomenclature of the corresponding complexes are given in the bottom. group. A coarse grid with 1 ˚ A grid spacing and a finer grid with 0.25 ˚ A grid spacing were used. Aspartates, glutamates, lysine, histidines (two sites for each histidines) cystein, tyrosine and Nand C-termini are treated as titratable groups. pKmodel avalues of the titratable groups are given in Table 6.1. The pKavalues calculated for the CuBcenter both in the reduced and in the oxidized states in the aqueous solution are taken as a model pKavalues of the CuBcenter in protein. 6.4.1 CALCULATION OF AVERAGE PKaIN PROTEIN The electrostatic calculations were performed on the protein structure of R. sphaeroides. The electrostatic energies obtained from the QMPB calculations were decomposed into Born, background and the interaction energies. The difference in the Born energy (∆∆GBorn) (see Eq. 2.36) and the background energy (∆∆Gback) (see Eq. 2.37) for a protonation reaction of 136 pKacalculations of CuBligands in cytochrome coxidase Titratable group pKmodel a aspartate 4.0 glutamate 4.4 arginine 10.4 lysine 10.4 histidine Nδ6.6 histidine N7.0 tyrosine 9.6 cysteine 9.5 C-terminus 3.8 N-terminus 7.5 Table 6.1. pKmodel avalues of the titratable groups. The pKmodel avalues are taken from Ref. [70] and Ref. [71]. the site can be directly obtained from the electrostatic calculations (see section 2.2.2). The average pKavalues of CuBligand in protein were obtained from Eq. (6.1). pKa,i =pKintr a,i + N X j=1 (hxji − x(0) j)Wij (6.1) where the intrinsic pKavalue (pKintr a,i ) is the pKavalue that the particular titratable group would have, if all other titratable groups are in their reference form. This term includes the solvation energy and the interaction with non-titrating residues and the protein backbone. The term Wij represents the interaction energy between the titratable groups iand jin their charged from; hxjirepresents the protonation probability of the group jwhich is obtained by a thermodynamic average over all possible protonation states and x0 jis the reference protonation form of site j. The hxjiare obtained from MC calculations. These average pKavalues do not represent an equilibrium situation, but they are a good approximation of the real pKa value of the site. The protonation probabilities were obtained by Metropolis Monte Carlo calculations. The titration curves were calculated by Monte Carlo program, GMCT developed in our group. The protonation probabilities were computed at pH 7. The temperature was set to 300 K. The double flip and triple flips were set to 2 and 3 pH units respectively. The number of full MC scans were set to 30,000 and 100 equilibrium MC scans were performed before the MC full scan. 6.5. Results and Discussion 137 6.5 RESULTS AND DISCUSSION The pKavalues calculated for the Cu-bound ligands in various CuBand heme a3complexes in cytochrome coxidase are discussed in this section. The pKavalues discussed in this section are average pKavalues of the deprotonation reactions at pH=7. The effect of the protein environment on the CuBligands were studied by combining the DFT and continuum electrostatic calculations. The pKavalues obtained for the CuBcenter both in oxidized and reduced state in aqueous solution are discussed in detail in chapter 5. PKaVALUES OF MICROSCOPIC STATES OF CuBLIGANDS IN CYTOCHROME cOXIDASE. The pKa values of all possible proton equilibrium reactions of oxidized (Cu+2 B) and reduced (Cu+1 B)CuB center in the presence of different heme a3states of cytochrome coxidase are given in Table 6.2 and 6.3 respectively. The pKavalue of H2O ligand was calculated both in oxidized and reduced states of the CuB center. There is no experimental pKavalue available for H2O ligand in the CuBcenter. The pKavalues of H2O ligand obtained in the aqueous solution are used as model pKavalues. Compared to the pKavalues of the H2O in aqueous solution, the protein environment shifts the pKavalues of the H2O both in oxidized and reduced state of CuBcenter to higher values. The pKavalues are very high for the reduced state of the CuBcenter compared to the oxidized state. The pKavalues are in the range of 51-60 which show that the deprotonation of the H2O in the reduced state is highly unfavorable. The pKavalues calculated in the FeII-OH state are obviously high due to -1 charge on the heme a3center. The pKavalues increase with the change in charge state of heme a3center from +1→0→-1. The high pKavalues for deprotonation of H2O show that it is protonated within the physiological pH. There is a possibility that His333 another ligand of CuBcenter, can serve as a proton loading site for subsequent pumping. The pKavalues of His333 ligand obtained in the aqueous solution is used as model pKavalues. There is no experimental pKavalue available for His333 ligand in the CuBcenter. The low dielectric protein environment shifts the pKavalues of the His333 both in oxidized and reduced state of CuBcenter to higher values compared to results in aqueous solution. The experimental pKavalue for the deprotonation of imidazole, which produces an anionic imidazole, is approximately 14 [193, 194]. In present calculations, the pKavalues for His333 deprotonation are in the range of 37-48 in the oxidized state and 44.553.8 in the reduced state. The high pKavalues in the protein show that His333 is likely to be protonated within the physiological pH both in the oxidized and in the reduced state of the CuBcenter. The protonated Nδ1of His333 is hydrogen bonded to Thr352. As a result, the deprotonation of His333 may be significantly hindered. The pKavalue of His334 obtained both in aqueous solution and protein (hpKa,iiprot) with different heme a3ligand are given in Table 6.2 and 6.3 respectively. The pKavalues of His334 ligand obtained in the aqueous solution are used as model pKavalues. There is no experimental pKavalue available for His334 ligand in the CuBcenter. The protein environment shifts most of the pKavalues of the His334 both in oxidized and reduced state of CuBcenter. The pKavalues are shifted to higher values in the oxidized state of the CuBcenter and are in the range of 27.5-39. The pKavalues are even shifted to higher values in the reduced state of the CuBcenter. The pKavalue of 33.5 is obtained in the case of oxo-ferryl state (FeIV=O) whereas 144 pKacalculations of CuBligands in cytochrome coxidase Subunit Residue hxiat pH=7 FeIII-H2O FeIII FeII FeII-H2O FeIII-OH FeIV=O FeII-OH [+1] [+1] [0] [0] [0] [0] [-1] GLU-533 0.002 0.001 0.002 0.002 0.002 0.002 0.002 HIS-534 0.411 0.409 0.436 0.437 0.426 0.434 0.461 ASP-536 0.001 0.003 0.003 0.002 0.003 0.004 0.004 GLU-539 0.001 0.001 0.001 0.002 0.002 0.002 0.003 TRP-540 1.000 1.000 1.000 1.000 1.000 1.000 1.000 GLU-548 0.004 0.003 0.003 0.003 0.003 0.003 0.003 HIS-549 0.107 0.118 0.123 0.119 0.122 0.117 0.128 GLU-552 0.002 0.003 0.004 0.004 0.003 0.004 0.003 LYS-556 0.997 0.998 0.999 0.997 0.998 0.997 0.999 ARG-557 0.999 1.000 0.999 0.999 0.999 0.999 1.000 GLU-558 0.003 0.003 0.003 0.003 0.002 0.002 0.002 ASP-559 0.002 0.002 0.002 0.001 0.002 0.001 0.002 TRP-560 1.000 1.000 1.000 1.000 1.000 1.000 1.000 HEA a-PROPA 0.000 0.000 0.000 0.000 0.000 0.000 0.000 HEA a-PROPD 0.000 0.000 0.000 0.000 0.000 0.000 0.000 HEA a3-PROPA 0.000 0.000 0.000 0.000 0.000 0.000 0.000 HEA a3-PROPD 0.000 0.000 0.000 0.000 0.000 0.000 0.000 II GLU-31 0.019 0.015 0.016 0.015 0.021 0.016 0.015 ARG-35 1.000 1.000 1.000 1.000 1.000 1.000 1.000 HIS-55 0.004 0.003 0.005 0.003 0.005 0.004 0.007 TRP-56 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ASP-58 0.000 0.000 0.000 0.001 0.001 0.000 0.005 TYR-78 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TRP-81 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ARG-82 1.000 1.000 1.000 1.000 1.000 1.000 1.000 HIS-84 0.000 0.000 0.000 0.000 0.000 0.000 0.000 GLU-85 0.000 0.000 0.000 0.000 0.000 0.000 0.000 LYS-86 0.924 0.921 0.923 0.925 0.920 0.925 0.924 ARG-87 1.000 1.000 1.000 1.000 1.000 1.000 1.000 LYS-89 0.366 0.355 0.363 0.362 0.356 0.354 0.371 ARG-93 1.000 1.000 1.000 1.000 1.000 1.000 1.000 HIS-96 0.045 0.046 0.048 0.051 0.048 0.053 0.052 continued on next page 6.5. Results and Discussion 145 Subunit Residue hxiat pH=7 FeIII-H2O FeIII FeII FeII-H2O FeIII-OH FeIV=O FeII-OH [+1] [+1] [0] [0] [0] [0] [-1] GLU-101 1.000 0.999 1.000 1.000 0.999 1.000 1.000 TRP-104 1.000 1.000 1.000 1.000 1.000 1.000 1.000 GLU-128 0.011 0.013 0.012 0.012 0.014 0.014 0.016 GLU-131 0.011 0.011 0.011 0.014 0.009 0.012 0.009 ASP-133 0.001 0.001 0.001 0.001 0.001 0.001 0.001 LYS-137 0.999 1.000 1.000 1.000 1.000 1.000 1.000 TYR-141 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TRP-143 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TYR-144 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TRP-145 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TYR-147 1.000 1.000 1.000 1.000 1.000 1.000 1.000 GLU-148 0.030 0.027 0.031 0.034 0.029 0.031 0.031 TYR-149 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ASP-151 0.006 0.008 0.006 0.009 0.009 0.009 0.008 GLU-152 0.005 0.005 0.006 0.005 0.004 0.007 0.006 GLU-153 0.010 0.010 0.009 0.008 0.010 0.010 0.012 GLU-157 0.080 0.080 0.080 0.086 0.075 0.078 0.081 TYR-159 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ASP-169 0.004 0.004 0.005 0.006 0.005 0.005 0.005 ARG-171 1.000 1.000 1.000 1.000 1.000 1.000 1.000 GLU-175 0.002 0.003 0.002 0.002 0.002 0.001 0.002 GLU-177 0.011 0.009 0.009 0.010 0.011 0.012 0.007 GLU-182 0.013 0.013 0.011 0.012 0.015 0.017 0.015 TYR-185 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ARG-187 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ASP-188 0.001 0.001 0.001 0.001 0.001 0.001 0.001 GLU-189 0.272 0.271 0.262 0.263 0.254 0.268 0.274 ASP-195 0.001 0.000 0.000 0.001 0.001 0.001 0.001 LYS-204 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ASP-214 0.000 0.000 0.000 0.000 0.000 0.000 0.000 TRP-219 1.000 1.000 1.000 1.000 1.000 1.000 1.000 LYS-227 0.984 0.984 0.996 0.997 0.997 0.997 1.000 continued on next page 146 pKacalculations of CuBligands in cytochrome coxidase Subunit Residue hxiat pH=7 FeIII-H2O FeIII FeII FeII-H2O FeIII-OH FeIV=O FeII-OH [+1] [+1] [0] [0] [0] [0] [-1] ASP-229 0.000 0.000 0.000 0.000 0.000 0.000 0.000 ARG-234 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TRP-239 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ARG-241 1.000 1.000 1.000 1.000 1.000 1.000 1.000 GLU-243 0.000 0.000 0.000 0.000 0.000 0.000 0.000 ARG-244 0.998 0.999 0.998 0.999 0.999 0.998 0.999 GLU-245 0.040 0.043 0.040 0.039 0.042 0.041 0.042 TYR-262 1.000 1.000 1.000 1.000 1.000 1.000 1.000 LYS-268 0.975 0.981 0.979 0.977 0.980 0.975 0.977 GLU-272 0.013 0.016 0.017 0.015 0.014 0.014 0.017 GLU-273 0.007 0.005 0.004 0.005 0.006 0.004 0.007 TYR-275 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TRP-278 1.000 1.000 1.000 1.000 1.000 1.000 1.000 GLU-280 0.016 0.012 0.015 0.014 0.015 0.015 0.014 ARG-283 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TYR-287 1.000 1.000 1.000 1.000 1.000 1.000 1.000 GLU-288 0.054 0.055 0.053 0.052 0.053 0.057 0.052 III HIS-3 0.245 0.241 0.236 0.233 0.232 0.241 0.239 LYS-5 0.997 0.998 0.998 0.998 0.998 0.997 0.997 HIS-7 0.000 0.000 0.000 0.000 0.000 0.000 0.000 ASP-8 0.000 0.000 0.000 0.000 0.000 0.000 0.000 TYR-9 1.000 1.000 1.000 1.000 1.000 1.000 1.000 HIS-10 0.001 0.001 0.001 0.001 0.002 0.002 0.002 TRP-17 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TRP-35 1.000 1.000 1.000 1.000 1.000 1.000 1.000 HIS-37 0.460 0.462 0.486 0.493 0.491 0.490 0.510 TRP-42 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TYR-53 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TRP-58 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TRP-59 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ASP-61 0.004 0.004 0.004 0.003 0.005 0.005 0.004 GLU-65 0.000 0.000 0.000 0.000 0.000 0.000 0.000 continued on next page 6.5. Results and Discussion 147 Subunit Residue hxiat pH=7 FeIII-H2O FeIII FeII FeII-H2O FeIII-OH FeIV=O FeII-OH [+1] [+1] [0] [0] [0] [0] [-1] GLU-68 0.046 0.045 0.042 0.040 0.041 0.043 0.042 ASP-70 0.837 0.836 0.847 0.848 0.848 0.845 0.845 HIS-71 0.000 0.000 0.000 0.000 0.000 0.000 0.000 ARG-76 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ARG-80 1.000 0.999 0.999 0.999 1.000 0.999 1.000 TRP-81 1.000 1.000 1.000 1.000 1.000 1.000 1.000 GLU-90 0.969 0.975 0.979 0.977 0.981 0.975 0.984 TRP-97 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TRP-99 1.000 1.000 1.000 1.000 1.000 1.000 1.000 LYS-103 1.000 1.000 1.000 1.000 1.000 1.000 1.000 HIS-104 0.000 0.000 0.000 0.000 0.000 0.000 0.000 TYR-107 1.000 1.000 1.000 1.000 1.000 1.000 1.000 GLU-112 0.018 0.018 0.019 0.017 0.017 0.021 0.020 ASP-117 0.020 0.021 0.025 0.025 0.022 0.024 0.031 GLU-123 0.057 0.060 0.057 0.059 0.055 0.060 0.057 ASP-129 0.000 0.001 0.000 0.001 0.000 0.001 0.000 TRP-131 1.000 1.000 1.000 1.000 1.000 1.000 1.000 HIS-132 0.013 0.012 0.016 0.015 0.014 0.013 0.017 CYS-143 1.000 1.000 1.000 1.000 1.000 1.000 1.000 CYS-146 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TRP-150 1.000 1.000 1.000 1.000 1.000 1.000 1.000 HIS-152 0.000 0.000 0.000 0.001 0.000 0.000 0.000 HIS-153 0.035 0.035 0.031 0.034 0.035 0.037 0.034 HIS-157 0.417 0.433 0.430 0.408 0.413 0.416 0.419 GLU-158 0.072 0.069 0.068 0.073 0.073 0.065 0.071 ARG-161 0.998 0.997 0.998 0.997 0.998 0.998 0.997 ARG-162 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ASP-163 0.000 0.001 0.000 0.001 0.001 0.001 0.000 TRP-166 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TYR-184 1.000 1.000 1.000 1.000 1.000 1.000 1.000 GLU-185 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TYR-186 1.000 1.000 1.000 1.000 1.000 1.000 1.000 continued on next page 148 pKacalculations of CuBligands in cytochrome coxidase Subunit Residue hxiat pH=7 FeIII-H2O FeIII FeII FeII-H2O FeIII-OH FeIV=O FeII-OH [+1] [+1] [0] [0] [0] [0] [-1] HIS-188 0.982 0.983 0.978 0.985 0.981 0.983 0.982 TYR-198 1.000 1.000 1.000 1.000 1.000 1.000 1.000 HIS-209 0.000 0.000 0.000 0.000 0.000 0.000 0.000 HIS-212 0.029 0.025 0.021 0.022 0.019 0.025 0.016 CYS-223 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ARG-226 0.321 0.332 0.342 0.250 0.256 0.254 0.361 ARG-229 0.997 0.996 0.997 0.997 0.997 0.997 0.996 HIS-231 0.001 0.001 0.001 0.001 0.000 0.001 0.001 GLU-235 0.002 0.002 0.002 0.001 0.001 0.002 0.001 LYS-236 0.290 0.275 0.275 0.299 0.298 0.302 0.282 HIS-237 0.521 0.512 0.505 0.562 0.548 0.561 0.512 GLU-241 0.140 0.138 0.150 0.156 0.156 0.150 0.148 TRP-245 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TYR-246 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TRP-247 1.000 1.000 1.000 1.000 1.000 1.000 1.000 HIS-248 0.000 0.000 0.000 0.000 0.000 0.000 0.000 ASP-251 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TRP-254 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TYR-262 1.000 1.000 1.000 1.000 1.000 1.000 1.000 TRP-264 1.000 1.000 1.000 1.000 1.000 1.000 1.000 IV HIS-11 0.316 0.318 0.320 0.320 0.317 0.316 0.319 ASP-17 0.000 0.000 0.000 0.000 0.000 0.000 0.000 GLU-22 0.002 0.002 0.002 0.003 0.001 0.002 0.003 LYS-23 0.993 0.996 0.994 0.995 0.995 0.994 0.995 ARG-30 0.999 1.000 0.999 0.999 0.999 0.999 1.000 TRP-34 1.000 1.000 1.000 1.000 1.000 1.000 1.000 The pKavalue of 46.4 and 46.2 are obtained for His334 in H0and H states respectively. The high pKavalues show that His334 remains protonated in these states. The next deprotonation is expected in connection with the formation of the O state and a pKaof 40.7 is obtained for His334 which is associated with proton pumping event. The high pKafor the His334 show that it is protonated in O state. The O state was proposed [1] to evolve in two different ways. The protonation of the hydroxyl group of the CuBcenter leads to O∼state and a pKavalue of 27.5 is obtained for the deprotonation of His334 and in next possible step (OH)apKavalue of 6.6. Conclusions 149 40.4 is obtained for His334. The pKavalues of 45.4, 32.7 and 45.9 are obtained for His334 in E, Epand R states respectively. All these high pKavalues show that His334 is protonated in O, E, Epand R states. In the present study, high pKavalues were obtained compared to the pKavalues reported by Stuchebrukhov and Fadda et al. This difference may be due to how the pKacalculations were done. For example the model pKavalues used for the CuBligands and the reliability of the charges used for the electrostatic calculations. In present study the average pKavalues of the CuBligands were calculated by Tanford-Roxby approximation [77]. Stuchebrukhov and Fadda et al. calculated the pKavalues by taking only the two instances like protonated and deprotonated states of the titratable groups and not all instances were considered by them. In our calculations, we considered all instances of the titratable groups, for example four instances were considered for the tyrosines and tryptophanes. An extensive bench mark studies were performed to obtain more accurate charges and solvation energies. In the two studies [1–3, 5, 6, 182], the crystal coordinates of the CuBcenter was not fully relaxed during the geometry optimization which can lead to unrealistic bond lengths. Fadda et al. only optimized the protonated structures, while the deprotonated was not optimized i.e., same geometry was used for both protonated and deprotonated states. The homogeneous dielectric model of the protein was used in the calculations. The charges of the CuBcenter used in the calculations were found without regard to the reaction field i.e., not in a self-consistent manner. The protonation state of the protein is not properly defined in their calculations where the charges of the titratable groups are critical for electrostatic calculations. The protonation depends on the redox states of the protein. This issues were completely ignored in Fadda et al. work. 6.6 CONCLUSIONS The DFT in combination with continuum electrostatic calculations were used to calculate the pKavalues of the CuBligands, H2O, His333 and His334 in various CuBand Fea3-CuB complexes to find the role of CuBligands in the reaction mechanism of cytochrome coxidase. The His334 which acts as ligand for the CuBcenter was proposed [1] to be involved in the reaction mechanism of cytochrome coxidase. In present study all the possible redox combination of CuB-heme a3center of cytochrome c oxidase are considered and all microscopic pKavalues of H2O, His333 and His334 were calculated for different states of the CuB-heme a3center. The pKavalues of His334 in F→O, E→Eptransitions of the catalytic cycle were studied. Our calculations are based on the detailed study of all the pKavalues of the ligands in the CuBcenter considering all redox states of CuBand heme a3center, which allows us to analyze the effect of the redox state of CuBheme a3center. According to the proposed scheme [3, 5, 6], His334 undergoes deprotonation whenever the chemical proton enters the binuclear center and convert the hydroxyl ligand to water molecule. This model was proposed based on the pKavalue of His334 where the pKa value of His334 was shifted to lower pKavalue around 5 pKaunits. In the present study, the pKavalues of CuBligands are significantly shifted to higher values in protein compared to aqueous solution. The pKavalues of CuBligands are in the range of 15-60. The high pKavalues of His334 show that His334 is protonated during all steps of the catalytic cycle. The pKavalues of His334 in protein increase significantly compared to the 150 pKacalculations of CuBligands in cytochrome coxidase aqueous solution when the CuBcenter is present in the oxidized states which is inconsistent with the pKavalues reported by the Stuchebrukhov et al. The protein environment shifts the pKavalues of all complexes to higher values, independently of the redox state of the metals. According to the pKavalues of His334 the proton pumping model as suggested by Stuchebrukhov [1] might not be possible. CHAPTER 7 CONCLUDING REMARKS AND OUTLOOK This thesis is focused on the pKacalculations of CuBand heme a3center of cytochrome c oxidase using density functional theory and continuum electrostatic models to understand the role of the CuB-bound histidines in reaction mechanism of cytochrome coxidase. The cytochrome coxidase energetically couples the electron transfer reactions associated with the reduction of oxygen to water, to pump protons across the membrane. Although a vast amount of structural and functional information of cytochrome coxidase is available from experimental and theoretical data, the actual step of coupling the redox reactions to the proton translocation is poorly understood. Recently Stuchebrukhov et al. [1–3, 5, 6] suggested that the deprotonation of the CuBligand His334 plays a central role in the proton pumping mechanism of cytochrome coxidase. According to this suggestion, His334 deprotonates at its Nδ1when the CuBcenter gets oxidized. From our results, His334 is protonated during all steps of the catalytic cycle in cytochrome c oxidase and concluding that a simple deprotonation of His334 in the protein environment is impossible in an equilibrium situation due to the high pKavalue that this group shows when it is bound to CuBcenter. A role of this residue in the mechanism of proton pumping might not be possible. The pKacalculations on CuB-bound His333 and His334 in aqueous solution (Chapter 5) lead to the following conclusions: The pKavalues for the first deprotonation reactions of His333 and His334 are too high (∼15) to allow their deprotonation in aqueous solution. The experimental pKavalue for the deprotonation of imidazole is approximately 14 [193, 194]. The double and triple deprotonation reactions are not energetically plausible in aqueous solution due to their high pKavalues. The pKavalues in the the range of 15.2-21 demonstrate that His334 and His334 are protonated at physiological pH. The pKavalues of CuBbound His333 and His334 in oxidized state do not show any significant change of the pKavalue due to their coordination to Cu ion in aqueous solution. The results obtained from reduced CuBcenter revels a small increase in the pKavalues of imidazoles. The pKavalues of His333 and His334 in the reduced state of CuBcenter are even higher to allow their deprotonation in aqueous solution. The results show that the solvation energies needed for pKacalculations should be calculated by including the electronic polarization due to the solvent explicitly i.e., in self-consistent manner to obtain proper trend in the pKavalues of the CuBcenter. 151 152 Concluding Remarks and Outlook The pKacalculations on CuB-bound His333 and His334 in cytochrome coxidase (Chapter 6) at pH=7 lead to the following conclusions: The pKavalues of CuBligands are significantly shifted to higher values in protein compared to aqueous solution. The pKavalues of CuBligands are in the range of 15-60. The CuB bound His334 is protonated during all steps of the catalytic cycle proving that the Fe and Cu ion oxidation states do not lower the pKavalues of CuBligands in the protein. The pKa values of CuBligands are significantly shifted to higher values dependent on the redox states of the CuBand heme a3center. The protein environment increases the pKavalues of all complexes to higher values independent of the redox state of the metals. The His333 and His334 are likely to be protonated at physiological pH. The double deprotonation and triple deprotonation reactions in protein are energetically unfavorable. The pKavalues of His333 show that this residue is likely to be protonated in the protein and an involvement of this residue as proton loading site in the reaction mechanism of cytochrome coxidase can therefore be ruled out. These results are inconsistent with the proposed role of His334 as a key element in the pumping mechanism. The proton pumping model as suggested by Stuchebrukhov [1] with the involvement of His334 might not be possible. A few issues could not be completely solved in the framework of this study and should be the topic of future research: In the present work, Tyr280 which is cross-linked with His284 is replaced by methyl group in the CuBmodel. This cross-linked His-Try residue can be represent by cross-linked imidazolephenol model and the pKavalues and redox potentials of the Tyr280 can be calculated to understand the role of cross-linked Try280 in the reaction mechanism of cytochrome coxidases. To understand the redox coupled protonation reactions in cytochrome coxidases, the midpoint potential of the redox centers CuA, heme a, heme a3and CuBcan be calculated using accurate theory which can be used to study redox behavior of the cytochrome coxidase. In the present work, minimal models are used for the CuBand heme a3centers. The interactions between the CuBand heme a3centers in the protein are treated as electrostatic interactions. Large models can be used for the binuclear center including the CuBcenter with the two histidine ligands and the cross-linked His-Tyr residues and heme a3center with substituted porphyrin. High level quantum chemical methods can be performed on this complete binuclear center, to study the influence of the heme a3redox states on CuBligands. The crystal structures of both mammalian and bacterial cytochrome coxidases do not contain all water molecules, that could be involved in the proton transfer pathways. Computer algorithms can be used to determine the likely positions and model internal water molecules in the protein cavities to understand the hydrogen bond networks in proton transfer pathways. We plan to use Vol program (developed in our group) to place water molecules and to find the hydrogen bond networks in cytochrome coxidase. The complex coupling of electron and proton transfer reactions in cytochrome coxidases must be elucidated in detail. In the future, we plan to use the DMC (Dynamic Monte Carlo) program to study the charge transfer and proton transfer rates in cytochrome coxidase. BIBLIOGRAPHY [1] D M Popovi´ c and A A Stuchebrukhov. Electrostatic study of the proton pumping mechanism in bovine heart cytochrome coxidase. J Am Chem Soc, 126:1858–1871, 2004. [2] D M Popovi´ c and A A Stuchebrukhov. Proton pumping mechanism and catalytic cycle of cytochrome coxidase: Coulomb pump model with kinetic gating. FEBS Lett, 556: 126–130, 2004. [3] J Quenneville, D M Popovi´ c, and A A Stuchebrukhov. Combined DFT and electrostatics study of the proton pumping mechanism in cytochrome coxidase. Biochim Biophys Acta, 1757:1035–1046, 2006. [4] A A Stuchebrukhov and D M Popovi´ c. Comment on “Acidity of a Cu-bound histidine in the binuclear center of cytochrome coxidase”. J Phys Chem, 110:17286–17287, 2006. [5] J Quenneville, D M Popovi´ c, and A A Stuchebrukhov. Redox-dependent pKaof CuB histidine ligand in cytochrome coxidase. J Phys Chem B, 108:18383–18389, 2004. [6] D M Popovi´ c, J Quenneville, and A A Stuchebrukhov. DFT/Electrostatic calculations of pKavalues in cytochrome coxidase. J Phys Chem B, 109:3616–3626, 2005. [7] O M H Richter and B Ludwig. Cytochrome coxidase - structure, function, and physiology of a redox-driven molecular machine. Rev Physiol Biochem Pharmacol, 147:47–74, 2003. [8] M Wikstr¨ om. Idetification of the electron transfers in cytochrome oxidase that are coupled to proton-pumping. Nature, 338:776–778, 1989. [9] P Mitchell. Coupling of phosphorylation to electron and hydrogen transfer by a chemiosmotic type of mechanism. Nature, 356:301–309, 1961. [10] D Voet and J G Voet. Biochemistry. John Wiley and Sons, Ltd, Chichester, 3rd edition, 2004. [11] B Ludwig and G Suhatz. A two-subunit cytochrome coxidase (cytochrome aa3) from Paracoccus denitrificans.Proc Natl Acad Sci USA, 77:196–200, 1980. [12] S Iwata, C Ostermeier, B Ludwig, and H Michel. Structure at 2.8 ˚ A resolution of cytochrome coxidase from Paracoccus denitrificans.Nature, 376:660–669, 1995. [13] C Ostermeier, S Iwata, B Ludwig, and H Michel. Fvfragment-mediated crystallization of the membrane protein bacterial cytochrome coxidase. Nat Struct Biol, 2:842–846, 1995. 153