Article Information

• PDF (348 kb)
• Full Text with Large Figures

PubMed

• Articles by António M. Baptista
• Articles by Cláudio M. Soares

Related Articles

• Reorganization and Conformational Changes in the Reduction o...

  Reorganization and Conformational Changes in the Reduction of Tetraheme Cytochromes Biophysical Journal, Volume 89, Issue 6, 1 December 2005, Pages 3919-3930 A. Sofia F. Oliveira, Vitor H. Teixeira, António M. Baptista and Cláudio M. Soares Abstract Molecular dynamics simulation (MD) constitutes an alternative to time-consuming experiments for studying conformational changes. We apply MD on a redox system where experimental information exists for the fully oxidized and fully reduced states: tetraheme cytochrome . Instead of doing one simulation for each state, we apply 10 4-ns replicas for both states, which provides robust statistics to characterize the redox changes. Besides these long simulations, we perform 120 short ones (50ps), where an equilibrated oxidized state is perturbed to a reduced state. This allows the application of a nonequilibrium method, the subtraction technique, which makes it possible to characterize the different timescales of conformational changes. Reduction induces conformational changes in the N-terminus and on the loops spanning residues 36–42 and 88–93, which correlate very well with experiments, demonstrating the applicability of this methodology. We also analyze the effect of reduction on hydrogen bonds, solvent accessible surface and bound water, the changes being found to involve the hemes and propionate groups. Redox-induced protonation is also investigated, by protonating the propionates D from hemes I and IV. Although this change in the former does not have major conformational consequences, it induces in the latter conformational changes beyond the ones obtained with reduction. Abstract | Full Text | PDF (430 kb)

• Modeling Electron Transfer Thermodynamics in Protein Complex...

  Modeling Electron Transfer Thermodynamics in Protein Complexes: Interaction between Two Cytochromes c3 Biophysical Journal, Volume 86, Issue 5, 1 May 2004, Pages 2773-2785 Vitor H. Teixeira, António M. Baptista and Cláudio M. Soares Abstract Redox protein complexes between type I and type II tetraheme cytochromes from Hildenborough are here analyzed using theoretical methodologies. Various complexes were generated using rigid-body docking techniques, and the two lowest energy complexes (1 and 2) were relaxed using molecular dynamics simulations with explicit solvent and subjected to further characterization. Complex 1 corresponds to an interaction between hemes I from both cytochromes . Complex 2 corresponds to an interaction between the heme IV from type I and the heme I from type II cytochrome . Binding free energy calculations using molecular mechanics, Poisson-Boltzmann, and surface accessibility methods show that complex 2 is more stable than complex 1. Thermodynamic calculations on complex 2 show that complex formation induces changes in the reduction potential of both cytochromes , but the changes are larger in the type I cytochrome (the largest one occurring on heme IV, of ∼80mV). These changes are sufficient to invert the global titration curves of both cytochromes, generating directionally in electron transfer from type I to type II cytochrome , a phenomenon of obvious thermodynamic origin and consequences, but also with kinetic implications. The existence of processes like this occurring at complex formation may constitute a natural design of efficient redox chains. Abstract | Full Text | PDF (589 kb)

• Calculated pH-Dependent Population and Protonation of Carbon...

  Calculated pH-Dependent Population and Protonation of Carbon-Monoxy-Myoglobin Conformers Biophysical Journal, Volume 80, Issue 3, 1 March 2001, Pages 1141-1150 Björn Rabenstein and Ernst-Walter Knapp Abstract X-ray structures of carbonmonoxymyoglobin (MbCO) are available for different pH values. We used conventional electrostatic continuum methods to calculate the titration behavior of MbCO in the pH range from 3 to 7. For our calculations, we considered five different x-ray structures determined at pH values of 4, 5, and 6. We developed a Monte Carlo method to sample protonation states and conformations at the same time so that we could calculate the population of the considered MbCO structures at different pH values and the titration behavior of MbCO for an ensemble of conformers. To increase the sampling efficiency, we introduced parallel tempering in our Monte Carlo method. The calculated population probabilities show, as expected, that the x-ray structures determined at pH 4 are most populated at low pH, whereas the x-ray structure determined at pH 6 is most populated at high pH, and the population of the x-ray structures determined at pH 5 possesses a maximum at intermediate pH. The calculated titration behavior is in better agreement with experimental results compared to calculations using only a single conformation. The most striking feature of pH-dependent conformational changes in MbCO—the rotation of His-64 out of the CO binding pocket—is reproduced by our calculations and is correlated with a protonation of His-64, as proposed earlier. Abstract | Full Text | PDF (257 kb)

• …more

Copyright © 1999 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 76, Issue 6, 2978-2998, 1 June 1999

doi:10.1016/S0006-3495(99)77452-7

Biophysical Theory and Modeling

Simulation of Electron-Proton Coupling with a Monte Carlo Method: Application to Cytochrome c₃ Using Continuum Electrostatics

António M. Baptista^,  , Paulo J. Martel and Cláudio M. Soares

Instituto de Tecnologia Química e Biológica, Universidade Nova de Lisboa, 2781-901 Oeiras, Portugal

Address reprint requests to Dr. António M. Baptista, Instituto de Tecnologia Química e Biológica, Universidade Nova de Lisboa, Apartado 127, 2781-901 Oeiras, Portugal. Tel.: 351-1-446-9613; Fax: 351-1-441-1277.

Redox processes are among the key functional aspects of living organisms, being involved in such important phenomena as respiration and photosynthesis. Many factors can affect or modulate redox processes, one of the more important being the pH of the medium: redox proteins often show a strong dependence of midpoint redox potentials (E^half) on pH and a corresponding dependence of midpoint pK_a values (pK^half) on potential, a phenomenon known as the redox-Bohr effect (Papa et al,Xavier, 1985). Given the fundamental electrostatic nature of the two types of processes, the existence of such a dependence is hardly surprising; instead, it would be strange if no coupling were found, considering the strength and long-range character of electrostatic interactions. Despite this nonspecific chemical-physical origin of the redox-Bohr effect, it is likely that evolution has “explored” this aspect and developed more specific and functionally relevant electron-proton couplings, which are thus expected to be present in many redox proteins.

The study of biomolecular systems by theoretical computational methods depends, to a great extent, on how closely one manages to model the actual in vivo or in vitro physical conditions of the system being investigated. Even when we limit ourselves to the study of the more easily modeled equilibrium conditions typical of in vitro systems, as is the case here, the nature of the models should be carefully considered. In particular, special attention must be given to the choice of external parameters of the system, the most important of which are the temperature and pressure and, in the case of electron-proton coupling, the electrostatic potential and pH. (In statistical mechanical terms, this corresponds to the choice of a proper ensemble.) Hence, our system of interest is a protein whose protonatable and redox sites display equilibrium fluctuations of state, and not an isolated protein with a particular redox and protonation state. If we regard the protein as “static” from the viewpoint of electron and proton binding, we will obtain a distorted picture of its behavior; in particular, some entropic factors will be lost (see below).

The level of treatment given to the external parameters (temperature, pressure, pH, electrostatic potential) is in general different for each computational method, and a complementary approach using several methodologies may be necessary. There are two major methodological routes to the study of equilibrium biomolecular processes that are essentially electrostatic in nature, such as redox or (de)protonation processes: molecular mechanics (MM) and continuum electrostatic (CE) methods.

The MM approach treats the protein as a flexible molecule surrounded by moving water molecules, the electrostatic interactions being just one of the various types of interactions between the atoms; conformational changes can then be studied by using molecular dynamics (MD) simulations. Because a single redox and protonation state of the protein has to be used, redox and protonation changes can be studied not with pure MD, but rather by using MM-based free energy methods (Zwanzig, 1954,Mezei and Beveridge, 1986,Beveridge and DiCapua, 1989). MM-based free energy calculations have been used with some success to compute redox potentials (Churg and Warshel, 1986,Cutler et al,Langen et al,Mark and van Gunsteren, 1994,Alden et al,Apostolakis et al,Muegge et al,Soares et al) and pK_a values (Warshel, 1981,Russell and Warshel, 1985,Warshel et al,Lee et al,Del Buono et al) in protein molecules; as usual, calculated values are relative to model compounds or to changes in the protein (mutations). The major problem with MM-based free energy calculations is their high computational cost, which led to the development of approximate methods, usually involving partial rigidification and/or some form of linear response hypothesis (Warshel and Russell, 1984,Muegge et al,Levy et al). Nevertheless, when many changes have to be considered simultaneously, such as when treating multiple proton equilibrium (a protein with s protonatable sites has 2^s different states), the MM approach is too demanding (less so if one uses a pairwise approximation; Del Buono et al) and cannot be used to make a complete study of the simultaneous protonation and redox equilibria in proteins. Although the MD approach can be extended to treat the effect of binding (Baptista et al), the resulting methods are still too demanding for the type of study we are interested in here.

The CE approach treats the protein and solvent as two distinct dielectric media, the protein atoms being represented as point charges in a static (average) molecular structure; a counterion atmosphere surrounds the protein, and the whole system is described by the Poisson-Boltzmann equation (Sharp and Honig, 1990,Honig and Nicholls, 1995). The free energies associated with charge modifications are then computed as the electrostatic energy change in this two-dielectric model. If the linear form of the Poisson-Boltzmann equation is used, the energy thus obtained is by definition pairwise decomposable. (Although potentials on the counterion region may exceed kT/e, the use of the nonlinear form poses nonadditivity problems that can only be solved by introducing additional approximations (Vorobjev et al); as usual in CE multiple binding problems, we will assume that the linear form gives reasonable potential values at protein atoms.) Thus, besides its modest computational requirements (compared with MD simulations), this approach directly produces a set of additive free energy contributions that makes it possible to address the simultaneous change of many charge states in a protein molecule. These states have to be weighted by using some additional method, the most general of which is perhaps a Monte Carlo (MC) method (Antosiewicz and Porschke, 1989,Beroza et al). The most evident problem of the CE approach is the use of a rigid protein structure: either we restrict the method to (hypothetical) rigid molecules (Baptista et al), or we adopt the use of a relatively high protein dielectric constant (4–20 or higher) (Gilson and Honig, 1986,Antosiewicz et al,Demchuk and Wade, 1996,Warshel and Papazyan, 1998) whose theoretical justification is not entirely clear. This moderately high dielectric response of the protein region seems to help compensate for the rigidity imposed on the structure and is useful in semirigid MM-based calculations (Warshel and Åqvist, 1991). Obviously, large conformational changes cannot be properly modeled with this approach, especially if they result from charge-induced reorganization (small conformational fluctuations are probably reasonably accounted for by the use of a relatively high dielectric constant for the protein); unfortunately, this is a problem that also afflicts the aforementioned partially rigid MM-based methods. Despite this conformational problem and the somewhat empirical parameterization associated with it, the CE treatment seems to capture many of the important features of electrostatic free energy changes. Many satisfactory CE calculations have been made of (relative) redox potentials (Bashford et al,Gunner and Honig, 1991,Soares et al) and pK_a values (Bashford and Karplus, 1990,Bashford and Gerwert, 1992,Bashford et al,Yang et al,Demchuk and Wade, 1996) in proteins, and these calculations seem to be the only feasible route to the study of electron-proton coupling including the full protonation equilibrium (Lancaster et al,Soares et al,Kannt et al,Martel et al), even if that implies sacrificing some configurational aspects (of both the protein and the solvent). Obviously, only equilibrium thermodynamic aspects of the coupling can be analyzed.

The CE studies of electron-proton coupling have so far investigated either the effect of the redox state on CE-based pK_a calculations (Lancaster et al,Soares et al,Kannt et al) or, conversely, the effect of the protonation state on CE calculations of redox potentials (Martel et al). Hence, although the energetics were always treated within the CE framework, the weighting of states was not done in an integrated way. The primary aim of this work is to extend the weighting of states to CE redox calculations and, more importantly, to integrate the simultaneous weighting of the states of protonatable and redox sites. By reflecting both the pH and electrostatic potential of the solution, this is the procedure that most realistically reproduces the actual conditions in most in vitro and some in vivo studies. This type of approach makes it possible to analyze couplings between sites, including electron-proton coupling, in a way superior to the usual “static” approach based on direct interactions. It also makes it possible to analyze binding fluctuations and their possible role in in vivo transfer processes.

From a thermodynamic standpoint there are two main reasons to include weighting of states when treating the equilibrium of ligand binding in general. The first is concerned with the inability of a single state to reproduce the average effect of the ligands; e.g., the energetic effect of a site whose mean occupation is 0.5 obviously cannot be reproduced by considering the site to be either empty or occupied. But another, more subtle reason does exist: like the deconstraining of any other extensive property (e.g., energy), the deconstraining of the number (and position) of ligands creates an additional source of entropic effects. This new entropic term does not arise from protein or solvent configurational changes that may occur upon binding, which also contribute to the overall entropy, but rather from the empty/occupied fluctuations of the binding sites. Formally, the new entropic term is identical to the entropy of a system whose energetic states correspond only to different occupation states of its sites (i.e., a system without different configurational energy states for the protein and solvent; see below), so that it may be called occupational entropy. The change in occupational entropy in some processes may be an important fraction of the total entropy, so it can be misleading to establish a necessary association between entropy and configuration, as in the entropy interpretation proposed for c-type cytochromes by Bertrand et al, in terms of solvent configurational changes. In particular, occupational entropy may contribute significantly to the E^half or pK^half of a site, because the change of the occupation state of that site can affect the occupation states of neighboring sites. If the equilibrium of redox sites is not treated explicitly, as is usually done in “static” CE redox calculations, their contribution to the occupational entropy is entirely lost. The secondary aim of the present work is to compute some of these occupational entropy contributions for the selected model system and infer from the general magnitudes whether they can play an important role in redox and protonation equilibria in general.

One system in which the redox-Bohr effect has been well characterized (mostly in terms of equilibrium thermodynamic properties) is the tetraheme cytochrome c₃ of Desulfovibrio spp. (Santos et al,Coletta et al,Salgueiro et al,Salgueiro et al,Turner et al,Turner et al,Park et al). These sulfate-reducing bacteria have the ability to use H₂ as the sole energy source when in the presence of sulfate (Badziong et al,Badziong and Thauer, 1978), and cytochrome c₃ is proposed to be involved in the sulfate reduction pathway, possibly by acting as a mediator between a periplasmic hydrogenase and transmembrane proteins that would pass electrons to the cytoplasm, where the reduction pathway takes place (Louro et al,Louro et al); the two electrons from H₂ oxidation are probably transferred simultaneously, eventually together with the two protons (Louro et al,Louro et al). The structure of cytochrome c₃ has been determined for several Desulfovibrio species (Higuchi et al,Matias et al,Matias et al,Czjzek et al,Morais et al) showing a generally conserved structure with four heme groups with bis-histidinyl axial coordination. The presence of several redox centers in such a small protein (13.5–15kDa), displaying the redox-Bohr effect, makes cytochrome c₃ an ideal model system for the theoretical study of electron-proton coupling. In this work we have chosen cytochrome c₃ from D. vulgaris Hildenborough (DvHc₃) as a case study. Experimental (Saraiva et al,Messias et al) and theoretical (Soares et al,Soares et al,Martel et al) studies of DvHc₃ made so far seem to indicate propionate D of heme I as the protonatable site more strongly involved in the redox-Bohr effect, although other sites (His⁶⁷, N-terminus, propionate D of heme IV) also seem to be involved in the effect; these characteristics seem to be generally preserved in other Desulfovibrio species (Martel et al).

In the present work we start by presenting the theoretical basis for the treatment of the simultaneous binding equilibrium of protonatable and redox sites, as well as its consequences for the general concept of coupling between different binding sites; the concept of occupational entropy is also discussed. We then present an application of the proposed methodology to DvHc₃, aimed at complementing studies using other approaches (Soares et al,Soares et al,Martel et al); the study also serves as an illustration of the potential of the method. We analyze the equilibrium of electron and proton binding in this protein, in particular the question of electron-proton coupling and the molecular basis of the redox-Bohr effect. A major result of the study is the correct prediction of the complete heme order in the reduction pathway of DvHc₃, which was not possible with previous “static” calculations. Site-site couplings are then analyzed in terms of “dynamical” equilibrium quantities (correlations), instead of the usual and more limited “static” analysis based on direct interaction (free) energies. The involvement of protonatable sites in the biological redox-Bohr effect is studied, and the results seem to corroborate previous work. The concerted transfer of electrons and protons is also discussed and analyzed in terms of binding fluctuations. Finally, the contribution of occupational entropies to E^half and pK^half values (which is necessarily absent from the usual “static” calculations) is examined, and its magnitude is found to be significant. This application also shows that, by establishing a more direct and realistic connection to experimental parameters, the new proposed method leads to results whose analysis is both much simpler (no subjective choices of binding states are necessary) and complete (all energetic and entropic factors are taken into account) than the one possible by the usual methods.

The binding state of a protein with s binding sites can be represented by a vector n=(n₁, n₂, …, n_s), where n_i=0 or 1, depending on whether site i is empty or occupied. The probability of occurrence of (or the fraction of molecules with) a binding state n, in conditions of thermodynamic equilibrium, is then given by

(1)

(2)

If we have a method for computing the ΔG(n) values, Eq. (2) can be used directly to perform a sampling of the states n, using, e.g., an MC method (Antosiewicz and Porschke, 1989,Beroza et al,Yang et al). To perform a full study of the relevant external parameters, we have to scan both the pH and potential of the system and simultaneously sample the protonable and redox sites, using a “two-way” MC scheme like the one described in the Methods section. The ΔG(n) values can in principle be obtained by considering a thermodynamic cycle involving model compounds (Warshel, 1981,Bashford and Karplus, 1990), as long as we can compute the free energy changes of the reactions 0→n in both the protein and the set of model compounds. As discussed in the Introduction, CE methods are the only feasible approach to such a large number of states, as in pK_a calculations; this methodology is discussed in the Methods section.

One way of measuring the coupling between sites is to analyze the magnitudes of their interaction (free) energies, because this gives a measure of how much a change in the binding state of one site perturbs the binding state of the other in a direct way. Thus the coupling between sites, and in particular of electron-proton couplings, is usually analyzed in terms of these direct interactions (Lancaster et al,Soares et al,Kannt et al).

The direct interaction is obviously an important contribution to the coupling between two sites (and probably the determinant one in many cases), but this “static” approach neglects the existence of indirect effects through other sites whose binding states may change in a concerted and complex manner (allostery is another, more familiar source of indirect effects). The tendency of two sites to have either identical or opposed binding states is determined by both the direct and indirect effects between them, and by looking only at the former we can lose important information, as the following simple example illustrates. Consider a three-site system in which site 1 has strong destabilizing interactions with sites 2 and 3, which in turn have a weaker destabilizing interaction between them. When the ligand activity is high enough to induce doubly occupied states, the simultaneous occupation of sites 2 and 3 is obviously the preferred one; hence, those two sites will tend to have the same binding state, even though their direct interaction is destabilizing. The existence of this particular doubly occupied state, which may be involved in the biological function of the system, originates from indirect effects and could not have been suspected from their destabilizing direct interaction alone; instead of looking at isolated direct interactions, one has to look at the set of all direct interactions (as well as the intrinsic affinities) and calculate the resulting equilibrium and correlations. This simple example illustrates the importance of including indirect effects when analyzing the effective coupling between two sites and suggests that a more appropriate measure for that coupling is the statistical correlation between their binding states. This correlation is given for a pair of sites (i, j) as a function of the variances and covariance of their binding states (Mood et al):

(3)

The analysis of coupling between sites based on statistical correlations is not limited to pairs of sites, but can also be applied to measure the coupling between groups of sites; in particular, it can measure the coupling between total protons and total electrons. This is simply an alternative way of looking at the traditional concept of linkage (Wyman, 1964) (the so-called linkage relation is a statement about two alternative ways of computing the covariance of the total numbers of different ligand types).

In general, this type of coupling analysis based on correlations, which is superior to the usual one based on direct interactions, is only possible when the full redox and protonation equilibria are considered, as in the method proposed here.

The statistics obtained from the sampling of binding states provide a route for computing various equilibrium thermodynamic properties, including the occupational entropies mentioned in the Introduction. The complete specification of the microstate of the protein system comprehends its binding state n and its configurational state Γ (the configurational state is here understood as the coordinates and momenta of the protein and surrounding solvent). The entropy of this system can be written as (Hill, 1960)

(4)

(5)

(6)

(7)

(8)

(9)

The contribution of the occupational entropy to the free energy change of a process is a measure of how much the fluctuations of the states n_i are altered by the change in question. In particular, the process of occupying a previously empty site will in general affect the occupational fluctuations of the remaining sites, meaning that the entropic contribution to the E^half or pK^half value of a given site is not necessarily configurational, as is often assumed (Bertrand et al). The standard free energy change of the reaction of binding to site i is given by

(10)

(11)

(12)

From a computational point of view, and as noted in the Introduction, the CE approach used here can only perform an approximated configurational averaging, especially for the protein conformation, where at best the averaging reflects small-scale fluctuations of the conformation used; this will necessarily affect the occupational entropies computed. The protein conformation in a binding system can, in principle, be explicitly treated along the line presented by Baptista et al, with consequent improvements in occupational entropy calculations, but such an approach is far too demanding for the full analysis intended in this work.

The structure of DvHc₃ used for the calculations was the one published by Matias et al, with PDB code 1cth, more precisely molecule A of the asymmetrical unit. Crystallographic water molecules were discarded. This structure was obtained under oxidizing conditions and may not be fully representative of more reduced states.

The choice of coordinates for polar and aromatic hydrogen atoms is an important aspect in CE calculations, because incorrect placement may lead to incorrect hydrogen-bond networks or even to atomic overlap for some protonation states (Alexov and Gunner, 1997). Hydrogen atoms were added in several stages. First, hydrogen atoms that are covalently nonambiguous and have a unique conformation were added using GROMOS 87 (van Gunsteren and Berendsen, 1987), considering explicit polar and aromatic hydrogens (Smith et al). Next, this structure was used to add hydrogen atoms that are covalently nonambiguous but have a variable conformation; this was done using the program WHATIF (Vriend, 1990), which implements an empirical potential for hydrogen bonding and optimization methodologies for hydrogen-bond networks (Hooft et al). Finally, with the structure thus obtained, this same procedure was used to add hydrogen atoms that are covalently ambiguous and correspond to the titrating protons of carboxylic groups. The positioning of the hydrogen atoms obtained in this last stage was subsequently checked by visual inspection of hydrogen-bonding patterns and corrected when found necessary; in particular, all four different positions in the carboxylic plane were considered, whereas WHATIF only considers the two extended ones.

After building of the DvHc₃ structure as described above, CE calculations were made using Donald Bashford’s software package MEAD (Bashford and Gerwert, 1992), the methodology of which is based on a thermodynamic cycle involving protein and model compounds and has been described elsewhere (Warshel, 1981,Bashford and Karplus, 1990). Although this software was originally developed for the treatment of protonatable sites, the treatment of oxidizable sites is formally identical, as discussed in the Theory section. MEAD assumes that a positive charge is added upon binding, so that the redox process has to be expressed in terms of oxidizable instead of reducible sites (see the Theory section). These oxidizable sites can then simply be included in the set of protonable sites given as input to MEAD, classified as cationic (thus, their “neutral” state is the reduced form).

The consideration of the thermodynamic cycle involving model compounds requires a specification of the structures of these compounds. In MEAD the model compounds of amino-acidic sites are obtained by “carving out” from the protein the N-formyl-N-methylamide derivative of the amino acid residue where the site is located. In the case of propionate sites the model compound was taken as the corresponding acetyl group. The model compound of the heme sites consists of the heme group, excluding the atoms used in the model compound of the propionates, plus the side chains of the axial histidines, and the C_β and S_γ atoms of the attached cysteines.

For each site i MEAD requires a pK_i^mod value, which originally corresponds to the pK_a of its model compound in solution. The pK^mod values for the usual amino-acidic protonable sites were taken from Nozaki and Tanford, 1967 and the value for propionate sites was taken as the pK_a of acetic acid (Weast, 1984); these are shown in Table 1. In the case of redox sites the pK_i^mod should correspond to the reduction potential of the model compounds used in the calculations, more precisely to −eE_i^mod/(2.3kT) (because they are oxidizable sites). In the present case, the model compound of the heme sites does not correspond to a real molecule whose potential can be determined or estimated easily. A reduction potential of −220mV has been measured for an octapeptide bis-histidinyl derivative of the heme of cytochrome c, at pH 7 (Harbury et al). Because our model compound has the heme more exposed to the solvent than this derivative, its more strongly charged form (the oxidized one) is in principle more stable; the propionate sites, although probably ionized at pH 7, are probably too strongly solvated to have a significant effect. Hence we expect a reduction potential slightly more negative than −220mV. In any case, because all redox sites are of the same type, with the same E_i^mod value, we can simply set this value equal to zero and consider all potentials as relative to this reference value, instead of the value of the standard hydrogen electrode. The standard reduction potential of our model compound can be determined a posteriori by comparison with experiment (see below).

The set of charges for the protonatable sites was taken or adapted from the GROMOS 87 force field (van Gunsteren and Berendsen, 1987) with polar and aromatic hydrogens (Smith et al). Because all titrating protons have been positioned in the structure (see above), we can describe the deprotonated/protonated change as the addition of an explicit hydrogen atom. In some cases this is highly desirable, such as when one of the oxygen atoms of a carboxylic group is establishing hydrogen bonds with other groups and the other oxygen atom is free; in this case the proton should be added to the free oxygen atom. However, there are cases in which such an approach represents no advantage, as in the case in which titrable sites are totally exposed and are not establishing any obvious interaction with other groups. In this case, an average description of the deprotonated/protonated change may be more adequate. For this reason, a choice of average or explicit protonation was made based on visual inspection for Asp, Glu, Tyr, and propionate sites. All other titrable sites had an average description of protonation, which corresponds to always (never) using the positions of the titrating protons for cationic (anionic) sites; in these cases the titration corresponds only to a change on the partial charges of permanent atoms. The atomic partial charges used for cationic and anionic sites are shown, respectively, in Table 2,Table 3. The atomic partial charges for the heme sites were obtained from quantum chemical calculations, as described elsewhere (Martel et al), and are shown in Table 4.

The dielectric boundary between the protein cavity and the solvent is computed by MEAD by rolling a spherical probe over the protein molecular surface. The protein atomic radii defining the molecular surface were taken to be half of the minimum energy distance between like-atom pairs in the GROMOS 87 force field (van Gunsteren and Berendsen, 1987). The solvent probe radius was 1.4Å, which should be a reasonable spherical approximation of the water molecule, and the ionic exclusion layer thickness was set at 2.0Å, which approximates the sizes of counterionic species in physiological conditions (Gilson and Honig, 1988). The dielectric constant used for the solvent region was 80, the approximate value for bulk water at room temperatures. The dielectric constant for the protein interior was chosen as 15, based on previous studies with this same protein (Soares et al,Martel et al), in consonance with suggestions by other workers (Antosiewicz et al,Demchuk and Wade, 1996). As noted in the Introduction, the use of moderately high values for the inner dielectric constant is supposed to reflect the effect of conformational fluctuations, even though the question remains on empirical grounds. The ionic strength used in the calculations was 0.1M, a value approaching the in vivo values, and the temperature was 300K.

MEAD uses a finite-difference technique to solve the linear Poisson-Boltzmann equation. A two-step focusing method (Gilson et al) was used: a first calculation using a (80Å)³ cube with a 1.0-Å lattice spacing, centered on the protein, followed by a second calculation using a (20Å)³ cube with a 0.25-Å spacing, centered on the titrable site; for the model compounds calculations, the sides of the cubes were, respectively, 60Å and 15Å.

As output MEAD produces a set of intrinsic pK_a values, {pK_i^int}, and a matrix of site-site interaction energies, {W_ij}. In the case of protonatable sites, pK_i^int stands for the intrinsic pK_a of site i when all other sites j≠i are in their “neutral” state (the reduced state for redox sites, because they are treated as cationic); for redox sites it corresponds to the intrinsic reduction potential E_i^int in the “neutral” environment, and more exactly to −eE_i^int/(2.3kT). The pairwise term W_ij is the site-site interaction between i and j, more exactly the difference between the electrostatic energies of the like pairs (0,0) and (1,1), and those of the cross-pairs (0,1) and (1,0). These two sets of values contain all of the information necessary to compute the populations of binding states of the protein (see below).

As discussed above, redox sites can be easily included in the formalism of protonatable sites. Hence, as can be shown from the formalism of the Theory section, the probability of state n can simply be written as a generalization of the expressions previously derived by other workers for the protonation case (Tanford and Kirkwood, 1957,Bashford and Karplus, 1990,Yang et al):

(13)

For the sampling of binding states we have implemented an MC method that follows the Metropolis criterion to accept or reject trial modifications of state (Metropolis et al). The trial modifications consist of flips of the state of binding sites, 0⇆1, which automatically ensures that the stochastic matrix of trial moves is symmetrical, as required by the Metropolis algorithm (Allen and Tildesley, 1987). To avoid “bottlenecks” in binding space, where the system may become trapped, double flips were attempted for pairs of strongly coupled sites (W_ij≥2 pH units), as done by Beroza et al. Instead of random selection of individual and paired sites, and to speed up calculations, trial flips were sequentially attempted for all individual sites and pairs; this also corresponds to a valid Markov chain, leading to a proper sampling of the binding states (Hastings, 1970).

In usual calculations for protonatable sites, several MC runs are made within the pH range of interest. In the present case a grid of pH versus E (electrostatic potential) has to be considered instead. A main calculation was made in which the pH was varied from −5 to 25 in steps of 0.2, and E was varied from −500mV to 500mV in steps of 10mV (relative to the E^mod value of the heme sites; see above); the temperature was 300K. Each MC run consisted of 40,000 steps, and each step was defined as a full cycle of trial flips over the lists of individual and paired sites. Because the errors in fluctuations and correlations are always larger than those of average values, longer runs were made to compute properties involving those quantities, for selected values of pH and E; fluctuations and occupational entropies were computed using 2×10⁵ steps, and site-site correlations using 10⁶ steps.

Because two thermodynamic parameters are being examined simultaneously, pH and electrostatic potential E, the titration behavior of the system is fully represented by titration surfaces instead of curves.

The total redox titration surface of c₃DvH, shown in Figure 1a, displays a marked dependence on pH. The titration curves (i.e., the lines parallel to the E axis) are shifted to more positive E values as the pH decreases. This type of shift is what should be expected on electrostatic grounds: a lower pH means a tendency of the protonatable sites to be occupied, increasing the overall charge of the protein and thus facilitating the binding of electrons, so that reduction can occur at higher E values. The resulting sideways shift on the titration surface is the signature of the redox-Bohr effect. A more careful examination of the surface shows that several shifts occur: two pronounced shifts, in the pH ranges 1–8 and 10–15, and a much smaller one around pH 17. As will be seen below, these coincide with the regions where large protonation changes occur in the protein. This simply reflects the nonspecific chemical physical origin of the redox-Bohr effect alluded in the Introduction: whenever protonable and redox sites coexist in a protein, a thermodynamic coupling between them is to be expected because of the strength and long-range nature of electrostatic interactions. This will be more evident for redox sites, because they will usually have several neighboring protonatable sites that can affect them; because there are usually much fewer redox sites than protonatable ones, the effect may not be significant for many of the protonatable sites (see below).

Display large version of this figure

Figure 1

Redox titration surfaces for the total protein (a) and heme I (b).

The titration surfaces for the individual hemes are very similar to the surface of total titration; one of them is shown in Figure 1b. Although the exact location and extension of the shifts are slightly different for different hemes, they also occur at the pH regions of maximum protonation changes.

It should be noted that even though three shifts are present on the surfaces of Fig. 1, only one of them is of biological interest. The optimum pH for the growth of D. vulgaris is 6.5 (Badziong and Thauer, 1978), and a redox-Bohr effect on DvHc₃ is known to occur in the pH range 5–8 (Turner et al). Hence, only the first shift, in the pH range 1–8, should be of biological interest.

It is interesting to see how well the present calculation reproduces the experimental results at physiological pH. A direct comparison with experimental individual titration curves is not possible, but we can make the comparison with a simple binding model (1 protonable+4 redox sites) by Turner et al; this model reproduces quite well the experimental data it was fitted to and can be used to obtain titration curves. A comparison between the titration curves of that model and our calculations, at pH 6.6, is shown in Figure 2ab. The E scale of plot b has been shifted to obtain an optimal superposition of the total titration curves of both plots; our E=0, which refers to the reduction potential of the heme model compound in solution (see Methods), is seen to correspond to a standard value (relative to the standard hydrogen electrode) of about −265mV. This value is slightly more negative than the −220mV measured for a similar model compound (Harbury et al), as expected from desolvation arguments (see Methods).

Display large version of this figure

Figure 2

Comparison of experimental (a, c) and calculated (b, d) redox data for DvHc₃, at pH 6.6. (a, b) Titration curves of individual hemes (labeled solid lines) and the total protein (dashed line). (c, d) Populations of the molecules with 0–4 (heme-reducing) electrons.

The calculated titration curves for individual hemes are ordered as the experimental ones, i.e., the present methodology predicts correctly the relative order of reduction of the four hemes of DvHc₃ at physiological pH. Nevertheless, some important differences exist, specially for hemes I and II. Besides being too far apart, the calculated titration curves of these two hemes show a much more gradual variation than those of plot a. The high steepness of the experimental curves reflects the positive cooperativity experimentally observed between hemes I and II (Turner et al). Given the repulsive nature of the electrostatic interactions between the sites, the most likely explanation for this behavior is the existence of an allosteric conformational change induced by the change in redox state of hemes I and/or II. The existence of reduction-induced conformational changes in DvHc₃ has recently been observed in a NMR structural study (Messias et al) and in a MD simulation study (Soares et al). Because the CE approach adopted here cannot account for conformational changes (except perhaps for small fluctuations; see the Introduction), it is not surprising that we fail to predict the positive cooperative behavior. The consequence of this failure is that, by missing the sudden and almost simultaneous reduction of the two hemes, we are assigning a too high probability for molecules with two reduced hemes, as can be seen in Figure 2cd; plot c refers again to the model by Turner et al (where the molecules with 0–4 oxidized hemes are called stages). It is remarkable that even while missing this effect, the present methodology is able to predict the correct order of reduction for all four hemes, suggesting that the basic energetics of the system is satisfactorily described by the CE model. Of course, the excessive stabilization of molecules with two electrons in our calculations shows DvHc₃ to be less efficient in the uptake of two simultaneous electrons than it actually is. The evolutionary advantage of developing the supposed conformational change would have been the usual one for developing positive cooperativity: the fine-tuning of function within a very narrow range of physiological conditions.

A more synthetic way of expressing the effect of pH upon the redox titration of DvHc₃, and the usual way of characterizing the redox-Bohr effect, is in terms of the pH dependence of the E^half values. Figure 3a shows this dependence for the four heme sites; the curves correspond to the intersection of the individual titration surfaces (like the one in Figure 1b) with the plane z=0.5. The largest shifts occur where large protonation changes take place, as expected from electrostatic effects and discussed above, apropos of the titration surfaces. This illustrates again the nonspecific nature of the redox-Bohr effect, which inevitably appears whenever redox and protonable sites coexist on the same molecule. The order of reduction of the four hemes in the pH range 5–8 is consistent with Fig. 2. The E^half curves cross each other several times within the whole pH range; in particular, the order of the curves at very low pH is precisely the opposite of that at very high pH. The large number of crossings in the plot indicates a close balance of effects, which, given the successful prediction of reduction order, seem to be properly captured by the present methodology.

Display large version of this figure

Figure 3

E^half values for individual hemes (a) and isoreduction curves corresponding to states M1–M3 (b), as a function of pH.

In Figure 3b we show the curves resulting from the intersection of the total titration curve of Figure 1a with the planes z=1, 2, and 3. These loci correspond to macroscopic states with an average number of (heme-r