Charge Parametrization of the DvH-c

Subscriber access provided by Milner Library | Illinois State University

Article 3

Charge Parametrization of the DvH-c Heme Group: Validation Using Constant-(pH,E) Molecular Dynamics Simulations João Henriques, Paulo Jorge Costa, Maria J. Calhorda, and Miguel Machuqueiro J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/jp3082134 • Publication Date (Web): 30 Nov 2012 Downloaded from http://pubs.acs.org on November 30, 2012

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Charge Parametrization of the Dv H-c 3 Heme Group: Validation Using Constant-(pH,E ) Molecular Dynamics Simulations. João Henriques,† Paulo Costa,‡ Maria José Calhorda,† and Miguel Machuqueiro∗,† Centro de Química e Bioquímica and Departamento de Química e Bioquímica, Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal, and Departamento de Química, CICECO and Secção Autónoma de Ciências da Saúde, Universidade de Aveiro, 3810-193 Aveiro, Portugal E-mail: [email protected]

Phone: +351-21-7500112. Fax: +351-21-7500088

∗

To whom correspondence should be addressed de Química e Bioquímica and Departamento de Química e Bioquímica, Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal ‡ Departamento de Química, CICECO and Secção Autónoma de Ciências da Saúde, Universidade de Aveiro, 3810193 Aveiro, Portugal † Centro

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract We studied the effect of using different heme group charge parametrization methods and schemes (Merz–Kollman, CHELPG and single- and multi-conformational RESP) on the quality of the results produced by the constant-(pH, E) MD method, applied to the redox titration of Desulfovibrio vulgaris Hildenborough cytochrome c3 . These new and more accurate charge sets enabled us to overcome the previously reported dependence of the methods performance on the dielectric constant, ε, assigned to the protein region. In particular, we found a systematic, clear shift of the E mod towards more negative values than those previously reported, agreeing with an electrostatics based reasoning. The simulations showed strong coupling between protonating/redox sites. We were also able to capture significant direct and, especially, indirect interactions between hemes, such as those mediated by histidine 67. Our results highlight the importance of having a good quantum description of the system before deriving atomic partial charges for classic force fields.

keywords: protonation, reduction, potential, RESP, Poisson-Boltzmann, Monte Carlo

Introduction The coupling of redox processes with protonation events is known to be essential for living organisms, playing a crucial role in respiration, photosynthesis and other energy transduction systems. This type of interaction was named as the redox-Bohr effect 1 . It can be described as a strong pH influence on heme reduction potentials, corresponding to a thermodynamic dependence between the reduction of redox sites and the ionization of protonable sites, due to their electrostatic interaction. This effect is of potential physiological significance in cytochromes c3 , since it occurs even at physiological pH values. Therefore, the redox-Bohr effect is thought to allow (thermodynamically) correlated electron and proton capture in the same molecule, which constitutes a physiologically relevant mechanism of transferring the electrons and protons generated by hydrogenase from molecular hydrogen oxidation. While hemes receive the generated electrons, protonatable sites in the protein are more susceptible to be protonated by consequent shifts in their pK a values 2–9 . 2 ACS Paragon Plus Environment

Page 2 of 35

Page 3 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Throughout the years a great effort has been put towards understanding the molecular basis, regulation and structure-function relationship of protein electron transfer. Even though this has been very well characterized experimentally (mainly thermodynamically) in tetraheme cytochromes c3 2–5,8 , many researchers have turned to computational methods in order to understand its molecular basis and, in particular, to identify key events in the process 10–16 . With electrostatic interactions being a key factor behind the thermodynamics of a redox center, it is not surprising that earlier studies mainly focused on continuum electrostatics methods (e.g. Poisson-Boltzmann, PB) coupled with Monte Carlo (MC) sampling of the protonation/reduction states 10,12–15 . Others relied on molecular mechanics/dynamics (MM/MD), where protein conformation changes are explicitly treated for a fixed ionization (driven by protonation or reduction) 11,15 . When taken individually, PB/MC and MM/MD methods seem quite incomplete, as PB/MC methods are well adapted to describe electrostatic changes on non-flexible molecules, while MM/MD methods deal essentially with conformational aspects of molecular species. However, if taken together these distinct methods exhibit somewhat complementary characteristics, which were used by several authors, giving birth to a few constant-pH MD methods and subsequent extensions/ implementations 16–31 . The theoretical aspects behind most constant-pH MD methods can in principle also be used to address redox processes. Recently, Machuqueiro and Baptista introduced the constant-(pH,E) MD method 16 , an extension to the stochastic titration method 27 . With this method it is possible to sample both the conformation and protonation/reduction states of a protein on a solution with a given pH and redox potential (E) value. The authors observed a dependence of the quality of their results on the dielectric constant, ε, assigned to the protein. This should only be expected for rigid-structure PB/MC studies, where the structural reorganization must be entirely captured by the dielectric constant. Non-structural events, such as charge redistribution upon change in heme redox state, were proposed as the main cause for the poor results at low ε. Most current force fields use a fixed-charge model, each atom being assigned a single value for the atomic charge that is not affected by the local electrostatic environment. However, the accumulation of charge resulting from the non-uniform electronic distribution over the atoms of the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

system will depend not only on the type of bonds, but also on its general environment, which may either exert a further increase on the asymmetry of the electronic distribution or attenuate it. The inclusion of electronic polarizability on force fields would, in principle, solve this problem, making them much less system-dependent. Yet, the introduction of polarizability into available force fields has been inhibited by the high computational expense associated with calculating the local electrostatic field. Therefore, accurate determination of the characteristic partial atomic charges of a given system is crucial for the performance of the force field used in MM/MD 32 . There are essentially two routes for determining the parameters of the potential energy function, either fitting them to the results of ab initio quantum calculations or to experimental data. Most force fields include parameters generated with both approaches, rendering them semi-empirical. Despite this, and regardless of the parametrization protocol policy behind each force field, partial atomic charges are, for practical reasons, mainly obtained from quantum mechanical computations. However, there is no universally accepted procedure, as the concept of partial atomic charge remains non-consensual 33 . Partial atomic charges, unlike the electron density, are not a quantum mechanical observable, i.e. they can not be calculated from first principles. Hence, all charge derivation methods are ultimately arbitrary. They can, however, be compared on the basis of the electrostatic performance of their partial atomic charges for a given molecular system. GROMOS force field charge parameters may be determined from different levels of quantum theory and basis sets, as long as they are well supported and lead to a good electrostatic representation of the system 34,35 . Machuqueiro and Baptista 16 adopted Oliveira et al. 36 partial atomic charges for the redox centers, but this charge set was unable to cope with the significant charge redistribution resulting from high heme-heme interactions at low ε 16 . A possible problem may lie in the fact that it was derived from the computed molecular electrostatic potential (ESP) of the X-ray structure of one of the four hemes (heme I) present in DvH-c3 . Also, owing to the structural nature of the heme group, some atoms (i.e., Fe) are buried deep inside the molecular surface and their point charges can be wrongly determined, resulting in conformation dependent charge sets. In order to overcome these limitations, the electronic detail of the system can be increased by


Page 4 of 35

Page 5 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


performing full or partial geometry optimizations at a given QM level of theory for each heme structure at both redox states. Crystallographic X-ray structures are available for specific redox states, and hemes undergo redox-dependent conformation changes. Finally, considerable attention should also be given to the different charge derivation methods, especially to the ones which allow a multi-conformational analysis. The most commonly used methods for deriving partial atomic charges from quantum mechanically computed molecular electrostatic potentials are Merz-Kollman (MK) 37,38 , CHELPG (CHarges from ELectrostatic Potentials using a Grid based method) 39 , and RESP (Restrained ElectroStatic Potential) fitting 40 . In both MK and CHELPG methods, partial atomic charges are fitted to reproduce the molecular ESP at a number of points around the molecule. The ESP is one of the most obvious properties to reproduce in order to get proper partial atomic charges to model short to long range moleculemolecule interactions 33,40,41 . Those two methods mainly differ in the choice of the points where the electrostatic potential is calculated. ESP based methods are usually less sensitive to the level of theory and basis sets used in the quantum calculations, than methods based on an analysis of the wave function. The obtained charges are much more dependent on the conformation of the molecule. The representative partial atomic charges for flexible molecules should hence be computed as average values over several molecular conformations 40–42 . Furthermore, one needs to be aware that buried atoms, located far away from the points at which the ESP is computed, may yield fuzzy results. The RESP method was developed to tackle the aforementioned problems with the MK and CHELPG methods. This is done by including a penalty function in the least-squares charge fitting procedure, in the form of restraints on non-hydrogen atomic charges to a target charge. The objective of the restraints is to hold down the ESP derived charges to a lower magnitude with only a minimal decrease in the quality of the fit 40 . Apart from the possibility of restraining partial atomic charges to a certain value, the RESP fitting process allows the use of forced symmetry, and multiconformational ESP derived charge fitting. Unrestrained ESP derived charges (Merz–Kollman)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

and the charges obtained from the population analysis methods, e.g. Mulliken charges 43 , perform poorly when compared with RESP derived charges, for the simulation of biomolecules 44 . Taking into account the work by Oliveira et al. 36 , and for comparison purposes, we also decided to explore all three methods directly over the X-ray structures of reduced and oxidized heme I of DvH-c3 . Desulfovibrio vulgaris Hildenborough cytochrome c3 is a small (Mr 13,000) globular and monomeric tetraheme protein present in the periplasm of the sulfate reducing bacteria Desulfovibrio vulgaris Hildenborough, consisting of 107 residues plus four hemes covalently bound to cysteines in the polypeptide chain, together with bis-histidinyl axial ligation 45 . DvH-c3 plays an important role in the periplasmatic metabolism of molecular hydrogen in sulfate-reducing bacteria, acting as a mediator between the periplasmatic hydrogenase (the electron donor) and the high molecular weight cytochrome (the electron acceptor) 45,46 . DvH-c3 hemes are covalently held together in close proximity, exhibiting strong coupling between themselves and nearby acid/base groups. The heme reduction potential is thus dependent on the oxidation state of the other three hemes (redox interaction potential) and on the pH (redox-Bohr effect). Together with a very low residue to heme ratio, we obtain a tight packed structure where the heme reduction potentials are very close (within a range of 80 mV or 1.3 pH units). This makes DvH-c3 an extremely demanding test case for reduction potential prediction, considering that 1.3 pH units is close to the constant-pH MD method pK a prediction error reported previously 30 .

Computational Details and Methods Parametrization of the Redox Centers All quantum mechanical calculations were performed with the program Gaussian 03 47 , using the B3LYP functional and LANL2TZ(f) or 6-31G(d) basis sets for Fe and all other atoms, respectively 48–54 . The effective core potential LANL2TZ(f) was used for d6 and d5 iron atoms. All charge parametrizations were performed on DvH-c3 heme(s) I or I-IV, depending on the nature of the charge derivation scheme, i.e. single- or multi-conformational. DvH-c3 heme models (see 6 ACS Paragon Plus Environment

Page 6 of 35

Page 7 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1 in Results) consisted of the respective heme c group, the axial histidines side chains and covalently bound cysteine side chains of molecule A of the X-ray structure of this cytochrome (Protein Data Bank entry 2CTH) 55,56 . 2CTH has two molecules in the asymmetric unit with almost identical conformations, therefore we only used one in this study. Both histidine and cysteine side chains are separated from the rest of the respective amino acid at the Cβ , which is treated as a methyl group in the model. Heme propionates were treated as independent titrating sites and, therefore, were not included in the heme model structure for charge parameterization. Each heme group is similar in chemical formulation, differing uniquely in its (X-ray) conformation. The open source chemistry toolbox OpenBabel was used to automatically introduce all hydrogens in each model 57 . Since cytochrome c3 heme-iron is low spin in both Fe(II) and Fe(III) states, all systems were treated as having a spin multiplicity of 1 and 2, respectively 58,59 , corresponding to an overall charge of zero and 1. Six different charge derivation schemes were performed. CHELPG, MK and single-conformational RESP fitting procedures were performed over DvH-c3 heme I X-ray structure model without geometry optimization. A multi-conformational RESP fitting was also performed over the models of DvH-c3 hemes I-IV using the unrelaxed X-ray structures (RESP MC X-ray). Finally, the two remaining charge sets were computed with a multi-conformational RESP calculation on the DFT geometry-optimized models of DvH-c3 hemes I-IV (RESP MC Free Opt) or optimized with constraints on critical atoms, i.e. the axial histidine and covalent cysteine Cβ atoms (RESP MC Fixed Opt). With this last charge set, we aim to model the expected heme group nonplanarity by freezing the atom positions identified as its key determinants 60 . Each geometry optimization was followed by a frequency calculation to ensure that all optimized geometries represented a local/global minimum of the respective potential energy landscape. The RESP fitting procedure was done using a modified version of the AmberTools 1.3 RESP program, i.e. a patched standalone version of the RESP program obtained from http://q4mdforcefieldtools.org 61,62 . The patched version was chosen since the standard RESP program handles only a limited number of ESP points, which proved to be insufficient for the desired level of accu-



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

racy. All charges related to the apolar hydrogen atoms which are omitted in the united-atom MD heme group model were summed into their covalent heavy atom partial charge. Notice that in GROMOS 53A6, the united-atom model of the histidine residue side chain comprises two more hydrogen atoms than that from GROMOS 43A1 34,35 . This results in four more heme group charged atoms for the GROMOS 53A6 force field.

Redox Titrations All MM/MD and PB/MC simulations of DvH-c3 were performed using the constant-(pH,E) MD method 16 . The six, previously mentioned, heme charge parametrizations were included in the original GROMOS 53A6 force field. An appropriate RESP X-ray parametrization was also added to the GROMOS 43A1 force field (the number of hydrogen atoms is different between 43A1 and G53A6), for a force field effect study. Fourteen protonable/reducible sites, including hemes I-IV, their respective propionates A and D, histidine 67, and the N-terminal were titrated at the pH value of 6.6 and reduction potentials of -140, -180, -220, -240, -260, -280, -300, -320 mV. Each solvent relaxation and full MM/MD cycles were 0.2 ps and 2 ps long, respectively. All conditions were simulated with three replicates, 15 ns each. A total of 144 (53A6) plus 24 (43A1) simulations, 15 ns long each, were performed, yielding a total ∼ 2.5 µs of simulation. The heme redox titrations were computed by averaging at each redox potential the occupancy states of all four sites over the final equilibrated segment, i.e. the last 10 ns of each run. The data for each heme was fit to a Hill equation in accordance with Machuqueiro and Baptista 16 . MEAD 2.2.0 and PETIT 1.5 were used for PB and MC calculations, respectively 12,63 . PB atomic charges and radii were taken from each modified GROMOS force field. Dielectric constants of 2 and 80 were used for protein and solvent, respectively. Temperature and ionic strength were set to 300 K and 0.1 mol/dm3 , respectively. The PB equation was solved using the finite difference procedure with grid spacings of 0.25 and 1.0 Å. Each MC calculation was performed with 105 iterations. 8 ACS Paragon Plus Environment

Page 8 of 35

Page 9 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


MM/MD simulations were performed using a modified version of GROMACS 3.2.1 64,65 . DvHc3 was placed in the center of a rhombic dodecahedric box with periodic boundary conditions, filled with 5879 water molecules from a simple point charge (SPC) water model 66 . The minimum distance between images of DvH-c3 was always higher than 20 Å. The equations of motion were integrated with a 2 fs time interval and bond constraining was achieved by applying the LINCS algorithm 67 . Non-bonded interactions were treated using a twin-range cutoff scheme with values of 8 and 14 Å. The neighbor list was updated every 5 integration time steps or, in other words, every 10 fs. Long range electrostatic interactions were treated with the generalized reaction field method, in which a dielectric constant of 54 was used 68 . The ionic strength was set at 0.1 mol/dm3 . Pressure and temperature were kept constant during the course of each simulation through the use of Berendsen’s weak coupling methods 69 . Berendsen’s thermostat was set to yield a constant bath temperature of 300 K with a coupling constant, τ, of 0.1 ps. Pressure was kept constant at 1 bar, while the coupling constant was set to 2 ps. Isothermal compressibility was set to 4.5x10−5 bar−1 . Before each MM/MD simulation, an energy minimization of each system (fully reduced and fully oxidized) was performed. All bonds were constrained by LINCS. Different MM/MD initiations were done for each DvH-c3 heme redox state, charge set, and respective replicate. The first initiation was 50 ps long and all protein atoms were restrained (1000 kJ.mol −1 .nm−1 ). Velocities were generated with different seeds for each of the three replicates. In the second initiation only Cα were restrained. All analyses were performed using the GROMACS package and in-house tools. Structural representations were done using PyMOL 0.99rc6 and Chemcraft 1.6, while graphics were generated using Gnuplot 4.2 patchlevel 4 70–72 . The calculations of correlation-corrected errors for averages were computed using standard methods based on the autocorrelation function of the property measured to determine the number of independent blocks in the simulations. 73



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Results and Discussion Charge determination: electronic population An electronic population analysis of all QM calculations for charge set derivation was conducted. As can be seen in Figure 1, in reduced systems the HOMO is a t2g orbital localized on the iron atom and, upon oxidation, the spin density was found mainly over the hemic iron atom (1.071). No significant spin contamination was found in any tested oxidized heme model. Spin contamination is negligible if the value of < S2 > differs from s(s + 1) by less than 10% 33 . These results indicate that the B3LYP functional, with the used basis set is appropriate to describe the ground state electronic structure of the heme model.

Charge determination: geometry optimization The heme group charges used by Machuqueiro and Baptista in the implementation of the constant(pH,E) MD method were those previously published by Oliveira et al. in 2005 16,36 . The charges were obtained by a RESP fitting procedure of the electrostatic potential from the X-ray structure of DvH-TpI-c3 heme number one (for short: heme I X-ray). The buried Fe charges were left unrestrained, which may lead to poor predictions. These charges were applied to all four heme groups without taking into account that each heme group has its own unique conformation. Furthermore, each heme undergoes a conformational change upon oxidation, which leads to eight different heme conformations. Also, because the X-ray heme group conformation is not necessarily the same as in solution, geometry optimizations can be used to relax the structures and eliminate intramolecular crystal strain. In this case, we performed two different approaches: a free optimization (no positional freezing) and a fixed optimization, where the axial histidines and covalent cysteines C-β atoms were fixed to their X-ray position. Table 1 shows the RMSD values between hemes before and after geometry optimizations. It is clear that heme III presents the conformation most different from heme I X-ray (0.67 Å). Even though hemes I and II appear to be relatively close (0.39 Å), assigning the same charge set to all 10 ACS Paragon Plus Environment

Page 10 of 35

Page 11 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1: HOMO representation of the reduced d6 Fe(II) in Heme I X-ray model. Orbitals were rendered using Chemcraft 71 .



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

hemes regardless of their conformational differences can indeed be a significant approximation, especially when these charges are derived from a single X-ray structure, which is by no means representative of any heme group present in a protein in solution, at any given redox state. Thus, we need to use a less conformationally dependent charge derivation method, as CHELPG, or a multi-conformational approach, like RESP. All optimized structures were also compared to the X-ray structure of heme I, primarily showing how much each heme group deviates from its original X-ray structure upon geometry optimization. It also shows how different each heme group is from the heme I X-ray structure, used primarily by Oliveira et al. in their heme group charge parametrization 36 . In general terms, all optimized heme groups seem to deviate in the same proportion from their original X-ray structure. Heme IV is the outlier here, as it relaxes significantly more than the others. The redox state seems to influence structure optimization, but to a lesser extent. Full and partial optimizations, even when only 4 atoms out of 99 are restricted to their X-ray positions, constitute indeed two different approaches, with markedly different RMSD values.

Charge determination analysis Table 2 comprises all the partial atomic charges obtained for each charge set. The charges of apolar hydrogens were added to their respective heavy atoms charge, according to the heme group model present in united-atom force fields. Certain empirical rules should be followed by the different charge groups. For example, the partial atomic charges of the heme-iron atom (ferrous and ferric) should be positive, while nitrogen atoms will probably have a negative partial atomic charge. Furthermore, oxidized iron should in principle present higher positive charge than its reduced counterpart. In the previously published charge set 36 , the reduced heme-iron charge was 48% higher than its oxidized counterpart. This big discrepancy is hard to interpret and probably had an important influence on the results. Here, with the exception of MK charges, all iron and nitrogen partial atomic charges reflect formal oxidation states and electronegativity. MK charges for NE21 and NE22 (axial coordinating 12 ACS Paragon Plus Environment

Page 12 of 35

Page 13 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Structure RMSD crossing table. Optimized heme structures against heme X-ray (Å).

X-ray I

Fixed Free

II

Fixed Free

Fixed Free

0

-

-

-

-

-

-

0.212

-

-

-

Reduced

0.320

-

-

-

0.332

-

-

-

0.386

0

-

-

Reduced

0.441

0.221

-

-

Oxidized

0.443

0.232

-

-

Reduced

0.467

0.321

-

-

Oxidized

0.503

0.331

-

-

0.670

-

0

-

0.719

-

0.241

-

Reduced Oxidized

0.710

-

0.240

-

Reduced

0.750

-

0.277

-

Oxidized

0.746

-

0.266

-

0.486

-

-

0

X-ray IV

IV X-ray

0.228

X-ray III

III X-ray

Reduced

Oxidized

Free

Heme

II X-ray

Oxidized

X-ray Fixed

I X-ray

Reduced

0.320

-

-

0.321

Oxidized

0.370

-

-

0.317

Reduced

0.333

-

-

0.448

Oxidized

0.370

-

-

0.457



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 35

Table 2: Atomic partial charge for the six different charge sets. Atoms ordered by number and label. Single conformational charge sets were calculated on heme I and are represented simply by the method name (CHELPG, Merz-Kollman and RESP). Multi-conformational approaches were calculated on the four hemes and are represented by the acronym MC with the respective geometry optimization method. Atom 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61

FE NA NB NC ND CHA HHA C1A C2A C3A C4A CMA CAA CHB HHB C1B C2B C3B C4B CMB CAB CBB CHC HHC C1C C2C C3C C4C CMC CAC CBC CHD HHD C1D C2D C3D C4D CMD CAD CBS1 SGS1 CBS2 SGS2 CB01 CG01 ND11 HD11 CD21 HD21 CE11 HE11 NE21 CB02 CG02 ND12 HD12 CD22 HD22 CE12 HE12 NE22

CHELPG Reduced Oxidized 0.958 1.152 -0.469 -0.488 -0.235 -0.279 -0.320 -0.343 -0.498 -0.533 -0.258 -0.265 0.159 0.171 0.186 0.206 -0.064 -0.041 0.002 0.015 0.160 0.186 -0.010 0.029 -0.012 0.025 -0.231 -0.233 0.147 0.158 0.042 0.083 0.015 0.006 -0.109 -0.047 0.139 0.151 -0.016 0.023 0.442 0.384 -0.114 -0.080 -0.330 -0.308 0.186 0.183 0.182 0.202 -0.011 -0.021 -0.192 -0.119 0.167 0.171 -0.017 0.026 0.458 0.398 -0.098 -0.063 -0.371 -0.346 0.243 0.239 0.218 0.240 -0.063 -0.047 -0.023 -0.003 0.213 0.235 -0.015 0.021 -0.019 0.019 0.138 0.131 -0.414 -0.350 0.133 0.119 -0.415 -0.346 -0.015 0.022 0.268 0.282 -0.467 -0.440 0.363 0.377 -0.116 -0.115 0.077 0.087 0.305 0.306 0.017 0.033 -0.424 -0.427 -0.004 0.032 0.227 0.243 -0.427 -0.407 0.358 0.374 -0.196 -0.184 0.140 0.145 0.145 0.165 0.083 0.090 -0.218 -0.244

Merz-Kollman Reduced Oxidized 0.674 0.663 -0.366 -0.360 -0.124 -0.176 -0.240 -0.252 -0.338 -0.344 -0.026 -0.061 0.124 0.145 -0.138 -0.091 0.173 0.172 0.100 0.095 0.030 0.085 -0.085 -0.037 -0.097 -0.050 -0.117 -0.154 0.135 0.154 -0.109 -0.012 0.159 0.103 0.110 0.186 -0.151 -0.119 -0.080 -0.029 0.222 0.159 -0.080 -0.046 -0.246 -0.236 0.189 0.188 -0.022 0.027 0.207 0.166 -0.157 -0.086 0.064 0.097 -0.074 -0.021 0.279 0.222 -0.070 -0.036 -0.209 -0.217 0.203 0.207 -0.005 0.054 0.104 0.088 0.137 0.147 -0.116 -0.078 -0.087 -0.040 -0.092 -0.046 0.086 0.079 -0.293 -0.228 0.098 0.085 -0.296 -0.229 -0.084 -0.040 0.434 0.423 -0.481 -0.445 0.373 0.387 -0.449 -0.441 0.232 0.239 0.102 0.083 0.125 0.146 -0.015 0.042 -0.063 -0.020 0.344 0.333 -0.401 -0.371 0.355 0.370 -0.440 -0.412 0.256 0.259 -0.020 -0.018 0.182 0.194 0.074 0.097

RESP Reduced Oxidized 0.491 0.686 -0.331 -0.355 -0.137 -0.221 -0.168 -0.231 -0.351 -0.379 -0.032 -0.071 0.117 0.139 -0.122 -0.073 0.174 0.169 0.088 0.091 0.029 0.077 -0.080 -0.034 -0.094 -0.047 -0.112 -0.145 0.135 0.152 -0.068 0.013 0.110 0.077 0.177 0.218 -0.181 -0.117 -0.072 -0.025 0.174 0.123 -0.068 -0.035 -0.208 -0.214 0.171 0.173 -0.037 0.023 0.173 0.133 -0.104 -0.030 0.013 0.052 -0.064 -0.013 0.253 0.195 -0.065 -0.032 -0.213 -0.216 0.200 0.203 0.076 0.110 0.038 0.042 0.197 0.190 -0.155 -0.099 -0.078 -0.034 -0.098 -0.050 0.087 0.082 -0.275 -0.212 0.096 0.083 -0.282 -0.214 -0.068 -0.024 0.373 0.362 -0.428 -0.393 0.361 0.374 -0.429 -0.399 0.236 0.236 0.018 0.021 0.168 0.178 0.069 0.055 -0.072 -0.028 0.373 0.362 -0.428 -0.393 0.361 0.374 -0.429 -0.399 0.236 0.236 0.018 0.021 0.168 0.178 0.069 0.055

RESP MC X-ray Reduced Oxidized 0.477 0.666 -0.198 -0.215 -0.209 -0.291 -0.248 -0.310 -0.266 -0.218 -0.059 -0.041 0.101 0.117 -0.073 -0.059 0.116 0.112 0.131 0.139 -0.097 -0.046 -0.083 -0.039 -0.083 -0.038 -0.039 -0.069 0.126 0.138 -0.130 -0.062 0.186 0.175 0.080 0.091 -0.067 0.046 -0.088 -0.046 0.147 0.105 -0.050 -0.022 -0.159 -0.252 0.138 0.180 -0.051 0.078 0.222 0.149 -0.110 -0.063 0.049 0.062 -0.082 -0.033 0.252 0.253 -0.080 -0.053 -0.159 -0.144 0.140 0.141 -0.011 -0.022 0.108 0.118 0.097 0.117 -0.016 -0.048 -0.080 -0.035 -0.079 -0.037 0.075 0.070 -0.244 -0.187 0.078 0.082 -0.258 -0.222 -0.048 -0.001 0.224 0.229 -0.349 -0.328 0.356 0.373 -0.224 -0.216 0.191 0.192 0.010 0.029 0.110 0.114 -0.020 -0.031 -0.052 -0.010 0.224 0.229 -0.349 -0.328 0.356 0.373 -0.224 -0.216 0.191 0.192 0.010 0.029 0.110 0.114 -0.020 -0.031


RESP MC Fixed Opt Reduced Oxidized 0.480 0.666 -0.134 -0.153 -0.249 -0.256 -0.198 -0.239 -0.102 -0.113 0.018 -0.004 0.095 0.109 -0.155 -0.110 0.119 0.113 0.123 0.148 -0.096 -0.091 -0.072 -0.037 -0.069 -0.028 -0.091 -0.073 0.122 0.130 -0.052 -0.035 0.151 0.125 -0.012 0.059 0.031 0.036 -0.067 -0.025 0.133 0.083 -0.021 0.009 -0.182 -0.195 0.132 0.143 -0.062 0.022 0.203 0.164 -0.045 -0.006 -0.037 -0.041 -0.077 -0.042 0.172 0.156 -0.038 -0.002 -0.090 -0.072 0.129 0.136 -0.137 -0.130 0.124 0.146 0.132 0.130 -0.167 -0.127 -0.068 -0.032 -0.075 -0.034 0.068 0.060 -0.225 -0.157 0.065 0.070 -0.238 -0.195 -0.032 0.012 0.203 0.199 -0.345 -0.293 0.355 0.365 -0.217 -0.202 0.175 0.179 0.011 0.012 0.109 0.118 -0.027 -0.040 -0.034 0.004 0.203 0.199 -0.345 -0.293 0.355 0.365 -0.217 -0.202 0.175 0.179 0.011 0.012 0.109 0.118 -0.027 -0.040

RESP MC Free Opt Reduced Oxidized 0.478 0.676 -0.094 -0.144 -0.187 -0.286 -0.151 -0.235 -0.117 -0.159 0.030 0.017 0.092 0.105 -0.178 -0.138 0.120 0.127 0.131 0.136 -0.133 -0.076 -0.073 -0.035 -0.068 -0.030 -0.056 -0.088 0.117 0.134 -0.109 -0.011 0.175 0.122 -0.004 0.070 -0.025 0.034 -0.070 -0.026 0.133 0.073 -0.024 0.010 -0.126 -0.155 0.132 0.136 -0.142 -0.025 0.253 0.205 -0.051 -0.014 -0.037 -0.006 -0.095 -0.059 0.164 0.148 -0.034 0.000 -0.101 -0.110 0.129 0.136 -0.134 -0.083 0.135 0.131 0.117 0.128 -0.154 -0.115 -0.071 -0.027 -0.072 -0.032 0.071 0.063 -0.232 -0.160 0.062 0.074 -0.233 -0.195 -0.035 0.007 0.212 0.209 -0.363 -0.318 0.358 0.370 -0.204 -0.195 0.172 0.175 0.044 0.040 0.100 0.109 -0.067 -0.053 -0.037 0.003 0.212 0.209 -0.363 -0.318 0.358 0.370 -0.204 -0.195 0.172 0.175 0.044 0.040 0.100 0.109 -0.067 -0.053

Page 15 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


nitrogen atoms, belonging to the axial histidine residues) were sometimes found to assume positive values. MK Fe(II) and Fe(III) charges were also not in the right order, but the net difference between them was insignificant (1.1%). This was not completely unexpected because ESP based charge derivation methods usually yield weak charge estimates for deeply buried atoms. To overcome this problem, we used Mulliken iron charges for the RESP fittings. Despite their limitations, Mulliken charges still represent the best estimates of deeply buried atoms like heme-iron. Since they do not depend on the conformation, all RESP charge sets (single and multi-conformational) have similar iron charges depending on the heme redox state. In general, CHELPG charges follow expected trends. However, CHELPG showed a tendency to create bigger dipoles between atoms of opposite charge signal. From Tables S1 and S2 it is evident that CHELPG charges diverge significantly from all the other charge sets. The CHELPG scheme sampling method differs from that of MK, and ultimately from RESP which is based on the MK scheme. Multi-conformational approaches also yield different results from those given by single-conformation methods. For example, within the RESP approaches taken here, the resulting charges tend to diverge as the degrees of freedom allowed in the geometry optimizations are increased. Interestingly, the RMSD values between redox states shown in Table S3 are almost constant throughout the different charge sets. We also noted that the RESP fitting procedure altered the MK charges quite significantly. Figure 2 and Figure S1 depict atoms with atomic charges varying the most throughout the six charge sets considered in this work. Atoms closer to the coordination center present the highest average variation. This is significant, as these atoms are closely involved in the redox event that occurs at each heme-iron coordination center.

Redox Titration Following Machuqueiro and Baptista ’s approach 16 , a model of five pairwise interacting sites (4 redox +1 protonable) were successfully fitted to the experimental data published by Turner et al. 4 . More recently, Paquete et al. re-analized this data and published new thermodynamic parameters 15 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 2: Color representation of the average partial atomic charge variation. Atoms with partial charge varying the most between charge sets are represented in red. From orange to white, the atomic partial charge variation decreases. that do not differ qualitatively from the previous ones 74 . In order to compare our new results direcly with published results using this methodology, we used the original data 4 as a reference throughout this work. Figure 3 shows the total redox titration curves obtained using different charge sets with GROMOS 53A6 force field in constant-(pH,E) MD simulations of DvH-c3 . In order to compare with previous work 16 , we also present the redox titration curve obtained with the GROMOS 43A1 force field, using a new RESP charge set calculated on heme I X-ray (43A1 RESP). As mentioned in the Computational Details and Methods section of this work, all redox titration curves were obtained using a protein dielectric constant of 2. Once again, following a reported methodology 16 , the reduction potential of the heme model compound (E mod ) was taken as the value obtained in a rigid-structure PB/MC study, i.e. -249 mV 14 . The E mod value corresponds to the reduction potential of a hypothetical heme model compound in solution. An experimental E hal f (approx. -220 mV) was determined by Harbury et al. in 1965 75 , for an octapeptide bis-histidinyl derivative of the cytochrome c heme group. However, this heme group is affected by the particular environment and can be regarded at best as an indicative value for our model compound. Teixeira et al. rigid16 ACS Paragon Plus Environment

Page 16 of 35

Page 17 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


structure PB/MC study showed that an E mod value of -249 mV yielded the lowest RMSD between the E hal f values obtained from experimental and theoretical methods 14 .

Figure 3: Redox titration curves for each charge set, using E mod = −249 mV 14 . The redox titration curve corresponding to the simulations at ε = 2 from Machuqueiro and Baptista 16 is also shown in a dashed black line for comparison of results. From Figure 3 one can clearly see the different behavior between past and present results, showing the sensitivity of the redox titration towards the different charge sets studied. Remarkably, we found a small horizontal displacement between the curves of all charge sets tested in this work. Compared to the redox titration curve extracted from Machuqueiro and Baptista 16 , our model seems to be titrating at a more plausible reduction potential range (assuming the E mod value of -249 mV 14 ). The proximity between our titration curves and the experimental model titration range should have a great impact on the fitted E mod and E hal f values 4 . These values, especially 17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the E mod , should agree better with each other and with experimental data, than previous published work. The charge set has a bigger or equivalent impact on the redox titration behavior than the force field used, as evidenced by comparing these results with those in the previous study using the same force field (GROMOS 43A1). Even though the partial charge derivation scheme is similar (both used a RESP fitting of the ESP of heme I X-ray structure), the quantum mechanical settings and parameters used in the computations were different. Another important difference was the use of the Mulliken partial charge for the buried iron atom. On the other hand, the redox titration curves obtained using two different force fields agree well with each other. The redox titration curves obtained for the simulations with GROMOS 43A1 (previous and present work) present a lower slope than those simulated with GROMOS 53A6. This becomes clearer when the total titration curves are fitted to the experimental model 4 , as shown in Figure 4. This was first reported in the previous work 16 , being attributed to excessive electrostatic interactions between hemes, a consequence of the low protein dielectric constant used, ε = 2. With the new charge sets, the excessive electrostatic interactions seem to be attenuated and the difference in the curves between force fields can also have a noticeable contribution from the intrinsic non-bonded/electrostatic differences in GROMOS 43A1 and 53A6 force fields 76 . Even though our results lack sampling at some regions of the titration curves, most charge sets lead to good fit to experimental data 4 . Among the GROMOS 53A6 simulations, the multiconformational RESP curves seems particularly good, the worst being relative to the CHELPG charge set. Table 3 comprises the E mod and individual heme midpoint reduction potentials (E hal f ) determined for each charge set and respective force field, and the RMSD of the computed E hal f values relative to the experimental model. Since there is no experimental E mod value for our heme model compound, we determined an optimal E mod value for each charge set tested according to Machuqueiro and Baptista 16 . This was done by shifting the E values and determining the RMSD between the four computed E hal f and those obtained experimentally, at each step. The E value which gives the lowest RMSD is taken as the best estimate of the optimal E mod .


Page 18 of 35

Page 19 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

(b)

(c)

(d)

(e)

(f)

Figure 4: Total titration curves of each charge set fitted to the experimental model 4 . (a) 43A1 RESP (b) RESP (c) RESP MC X-ray (d) RESP MC Fixed Opt. (e) RESP MC Free Opt. (f) Merz-Kollman (g) CHELPG. 19 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 35

(g)

Figure 4: Total titration curves of each charge set fitted to the experimental model 4 . (a) 43A1 RESP (b) RESP (c) RESP MC X-ray (d) RESP MC Fixed Opt. (e) RESP MC Free Opt. (f) Merz-Kollman (g) CHELPG.

Table 3: RMSD values and corresponding E mod and E hal f obtained for each charge set. Experimental values are presented as a reference. Machuqueiro and Baptista’s 2009 results are relative to their simulations performed at ε = 2. For the Hill coefficient values, see Supporting Information. E hal f (mV)

RMSD

E mod (mV)

heme I

heme II

heme III

heme IV

range

mV

pH units

-126

-360

-316

-243

-283

117

57

0.96

43A1 RESP

-282

-252

-369

-219

-360

150

88

1.48

RESP

-304

-296

-320

-251

-334

83

58

0.98

RESP MC X-ray

-327

-303

-325

-229

-345

116

70

1.18

RESP MC Fixed

-320

-312

-304

-236

-350

114

69

1.16

RESP MC Free

-322

-306

-318

-242

-334

92

61

1.03

MK

-304

-306

-311

-235

-351

116

69

1.16

CHELPG

-294

-330

-297

-245

-329

85

60

1.01

Experimental 4

-220 75

-302

-307

-336

-256

80

Charge set Machuqueiro and Baptista

16


Page 21 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The computed E mod values agree well with one another, varying smoothly around ≈ -300 mV. However, they disagree greatly with the previously determined value, i.e. -126 mV 16 . A shift of the E mod towards more negative values than the E hal f obtained by Harbury et al. in their system (-220 mV) 75 is highly desirable. The reason is that our model compound heme groups are more exposed to the solvent than Harbury’s octapeptide model. Considering that the oxidized form of our model has a positive global charge, it should in principle be more stable in solution than the respective octapeptide bis-histidinyl derivative counterpart. Furthermore, the propionate sites, although probably ionized at the simulated pH, are very likely to be strongly solvated, having no significant effect on the E mod 12 . As a consequence, we expected the reduction potential of our model compound to be significantly more negative than the reported value, -220 mV 12,75 . As for the E hal f values, we noted that the values relative to each heme for each charge set are in close agreement with one another. Hemes I and II seem to titrate at relatively close reduction potentials while hemes III and IV show strong divergence in their E hal f values. This is reflected on the heme E hal f range. With the single exception of 43A1 RESP, the absolute E hal f range values lie in between 83 and 116 mV. These values are in close agreement with the previously reported range for the heme reduction potential, simulated at ε = 2 (see Table 3 first row) 16 . However, they assume marginally higher values than that reported by Turner et al. (see Table 3 last row) 4 . From the RMSD, we observed that each simulated charge set (with the exception of 43A1 RESP) yields similar results. The stochastic titration method seems to be rather robust to small changes in heme group partial atomic charge parameters. Moreover, the computed RMSD values, when translated into pH units, are close to that of the published stochastic titration method precision error, ∼ 0.8 pH units 30 . This emphasizes the high accuracy of the overall methodology employed in these new studies. Table 4 presents the order of reduction of the heme groups for each tested charge set, obtained by sorting the E hal f values of Table 3, in ascending order of E. As can be seen, previous results were not able to predict the correct order 16 . With the exception of RESP MC Fixed, CHELPG and 43A1 RESP we are able to predict correctly the relative order of reduction of heme I and II.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 35

Hemes III and IV are in the wrong order, as also happened in previous work by Machuqueiro and Baptista , even using higher dielectric constants 16 . Indeed, only when simulating at ε = 15, were the authors able to obtain the correct order of reduction for heme IV relative to heme III. With the exception of CHELPG and 43A1 RESP, which overestimate the negative nature of the reduction potential of hemes I and II, respectively, our results show that heme IV systematically titrated at lower E and heme III at higher E. A clear pattern can be extracted from Table 4 as E increases: IV < II ≤ I < III. The systematic nature of these results is remarkable, but, regardless of the charge set being used, calculated and experimental data are in complete disagreement when predicting the E hal f values for hemes III and IV. Even though we are aware of the difficulty to assign the E hal f values using NMR spectroscopy 5 , we found no ambiguity in this assignment for hemes III and IV. Table 4: Calculated and experimental order of reduction of the heme groups for each charge set. Each correct heme reduction order is shaded in gray. Charge set Machuqueiro and

Baptista a

43A1 RESP

low E

high E

I

Charge Parametrization of the DvH-c

Recommend Documents