INTRODUCTION

Top Abstract Introduction Methods Results & Discussion Conclusions References

The molecular mechanisms behind redox processes in proteins are largely unknown. There are questions that may be placed at several levels, namely the electron transfer process itself through the protein, the mechanism by which the protein controls the redox potential of the redox center, and the molecular motions due to the redox process. In theory, molecular dynamics simulation techniques may be able to investigate the last two questions by being capable of studying the dynamics and thermodynamics of diverse processes in proteins and other molecules. In practice, however, the consequences of a certain number of approximations that have to be introduced in the methodology are not obvious. For instance, electrostatic interactions may be one of the main factors determining the thermodynamics of a redox center and it is therefore important to know whether these interactions are sufficiently well implemented in a common force field to capture the reality of the redox process.

Despite the fact that continuum electrostatic (Gunner et al., 1997; Gunner and Honig, 1991; Lancaster et al., 1996) and protein dipoles Langevin dipoles (PDLD) (Alden et al., 1995; Churg and Warshel, 1986; Jensen et al., 1994; Warshel et al., 1997) methods are much more common in the study of this kind of systems, the use of molecular mechanics/dynamics techniques is a way to complement and possibly go beyond these static and semistatic approaches. There are some examples of studies where different redox states of protein systems have been analyzed using molecular mechanics/dynamics procedures (Leenders et al., 1994; Shenoy and Ichiye, 1993; Yelle et al., 1995). Similarly, there have been studies aiming to quantify thermodynamic quantities, such as redox potentials or redox potential differences, in organic molecules (Compton et al., 1989; Reynolds, 1990; Reynolds et al., 1988; Wheeler, 1994) and in proteins (Alden et al., 1995; Apostolakis et al., 1996; Mark and van Gunsteren, 1994).

Cytochrome c₃ is a small tetraheme protein present in the periplasm of Desulfovibrio species (see Coutinho and Xavier, 1994 for a review). The structure of this cytochrome has been solved (Czjzek et al., 1994; Higuchi et al., 1984; Matias et al., 1993, 1996; Morais et al., 1995) for several species. All four hemes are covalently bond to the protein and have bis-histidinyl axial coordination (Fig. 1). It has been shown (Badziong and Thauer, 1978; Badziong et al., 1978) that in the presence of sulfate, Desulfovibrio can use H₂ as its only source of energy. Under these conditions it has been proposed (Louro et al., 1997; Xavier, 1986) that cytochrome c₃ acts as a coupling protein to periplasmic hydrogenase being able to receive the two protons and two electrons resulting from the cleavage of H₂ (Louro et al., 1996, 1997). This concerted transfer of electrons and protons has been named the redox-Bohr effect (Papa et al., 1979; Xavier, 1985). Cytochrome c₃ is a quite interesting model system for studies of redox processes in proteins for several reasons: despite its small size, cytochrome c₃ displays complex redox behavior due to the presence of four closely interacting hemes in the same molecule, and due to the fact that each heme is in a different microenvironment. For cytochrome c₃ from Desulfovibrio vulgaris Hildenborough, whose structure (Matias et al., 1993) is used in this work, there is a large amount of experimental data, especially on the thermodynamics of reduction of all possible redox states (Salgueiro et al., 1992; Turner et al., 1994, 1996).

View larger version (49K): [in this window] [in a new window]   FIGURE 1   Ribbon representation (Kraulis, 1991) of cytochrome c3. The four covalently bound hemes presented in the molecule are depicted with balls and sticks.

The objective of this work is to use molecular dynamics simulation to study the effects of reduction upon structure dynamics. Another goal of this work is to apply and develop free energy calculations of the relative reduction energy of each individual heme. To achieve this goal, specific techniques to improve the accuracy of electrostatics in free energy calculations have been applied and tested, showing a large improvement in the results.

	MATERIAL AND METHODS

Top Abstract Introduction Methods Results & Discussion Conclusions References

General setup of molecular dynamics simulation

The molecular dynamics simulations described in this paper were performed with the package GROMOS 87 (van Gunsteren and Berendsen, 1987) with some modifications in order to implement some methods described below. The force field parameters used were not the original set of GROMOS 87 parameters, but a modified version (Smith et al., 1995) where aromatic hydrogens are explicitly treated and containing modified van der Waals parameters for the interaction of water oxygens and CH_n groups. The parametrization of heme groups had to be changed from the original parameters, which were meant for a hemoglobin type heme, while cytochrome c₃ contains c-type hemes covalently linked to cysteines. The vinyl groups with aromatic character (CAB-CBB and CAC-CBC) had to be changed to aliphatic groups CH1CH3 where the attachment of the cysteines takes place (CAB and CAC). To do this, changes in bonds, bond angles, and dihedrals had to be made in order to obtain the proper geometries. In keeping with the general philosophy of the force field, four aromatic hydrogens were added to the heme topology. The two cysteines and the axial histidines were then attached to the heme. The partial charges used on the heme atoms are in Table 1, where the changes between the reduced form (the form parametrized in the original force field) and the oxidized form are also listed. The covalently bond cysteines had the charges of C and S set to zero. The charges of the neutral histidine (Form A) of GROMOS were assigned to the axial histidines.

The cytochrome c₃ structure of Desulfovibrio vulgaris Hildenborough at 1.9 Å resolution (Matias et al., 1993) was used in the studies described here. Polar and aromatic hydrogens were added to the protein and to a set of chosen internal crystallographic waters. The water molecules included were establishing at least one hydrogen bond with atoms of the protein or they were completely buried and conserved in similar and homologous structures of cytochrome c₃ (Matias et al., 1993). According to pK_a calculations done on this system (Soares et al., 1997) and considering that the simulation is going to be done at pH 7, it was decided that bases are protonated, acids are deprotonated, the N-terminal is a NH₃⁺ group, and the only free histidine (His-67) is a neutral histidine deprotonated at N1. The position of the hydrogen atoms of the water molecules was optimized by energy minimization while keeping the protein rigid. This system was then solvated using a truncated octahedron with initial size of 58.884 × 58.884 × 58.884 Å³, generated by duplication of an initial equilibrated box containing 216 SPC water (Hermans et al., 1984) molecules. A minimum distance of 7 Å between the protein and box walls was imposed to generate the solvated system. The final system had 2792 water molecules and a total of 9615 atoms. This system was then energy minimized to remove excessive strain. This consisted of 1000 steps of steepest descent minimization of the water molecules followed by 1000 more steepest descent steps of the water and protein side chains.

The molecular dynamics simulations were performed using heat and pressure baths (Berendsen et al., 1984) at the constant temperature of 293 K, using a coupling constant of 0.1 ps for protein and water (independently), and at the constant pressure of 1 bar, with a coupling constant of 0.5 ps. For some of the simulations a special initialization procedure was used (see below), but in all cases the initial velocities were taken from a Maxwell-Boltzmann distribution at 293 K. SHAKE (Ryckaert et al., 1977) was used for all bonds with a geometric tolerance of 0.001. The integration of the equations of motion was carried out using a 0.002 ps time step. The twin range method for cutoffs (van Gunsteren and Berendsen, 1990) was applied to the charged-group pair electrostatic interactions: for charge groups at a distance below 8 Å the electrostatic force was calculated at every integration step, while for charge groups at a distance between 8 Å and 10 Å this interaction was only calculated each 10 steps. The heme redox charge groups were treated with a modified version of the twin range method, which is explained below.

Conformations along the molecular dynamics simulation were saved each picosecond. Average structures were obtained using conformations fitted to the first structure of the selected group of conformations, using C atoms. Hydrogen bonds were selected using a criterion of a maximum distance of 2.5 Å and a minimum angle of 135°. The NICE package (Nilsson, 1990) and our own programs were used in all the analysis presented here.

Calculation of the relative free energy of reduction

For calculation of the relative reduction free energy the well-known free energy perturbation techniques were used (Zwanzig, 1954). This particular type of problems has been analyzed before in other systems using these methodologies (Alden et al., 1995; Apostolakis et al., 1996; Compton et al., 1989; Mark and van Gunsteren, 1994; Reynolds, 1990; Reynolds et al., 1988; Wheeler, 1994). For the problem we are studying here, i.e., the calculation of relative free energies of reduction of individual hemes, the thermodynamic cycle represented in Fig. 2 a can be considered.

View larger version (19K): [in this window] [in a new window]   FIGURE 2   (a) Thermodynamic cycle that expresses the reduction process with two individual hemes. H1 and H2 represent two heme groups species. (b) Thermodynamic cycle that expresses the pseudoreduction process of one heme H1. The subscripts "g" and "prot" denote the gas phase and the protein, respectively; the superscripts "ox" and "red" denote the oxidized and the reduced species, respectively; F denotes the Faraday constant; n the number of electrons involved in the process, and E denotes the standard electrode potential.

Vertical energies in the cycle correspond to transfer energies of different species from the gas phase to the protein, and horizontal energies correspond to the complete reaction in the gas phase (G_g) and in the protein (G). The latter is proportional to the difference between redox potentials of the hemes considered. The free energy of the reaction in the protein can be written as:

(1)

Therefore, the free energy of the reaction (and the difference between redox potentials) can be calculated from the transfer energies of all species, plus the gas phase free energy difference. This is the general method for a calculation of this type, but the characteristics of the system studied here allow a significant simplification: if the redox system (the actual group that changes quantum mechanically) is confined at the heme and axial ligands, the gas phase free energy difference of the reaction is zero, given that all hemes in cytochrome c₃ are of the same type and have the same axial coordination. In other words, the species "transferred" to the gas phase are identical. This would be true for any calculation dealing with bihistidinyl coordinated c-type hemes, even if they belong to different proteins.

Therefore, only the transfer energies are required for calculating G. These can be calculated from individual pseudoreduction processes, such as the one described in Fig. 2 b, and the free energy relationship derived from it written in Eq. 2:

(2)

By using free energy perturbation techniques it is possible to simulate the (unphysical) process represented in the lower part of the cycle

(3)

and calculate its free energy change. The calculated free energy from the simulation is equal to the transfer energy terms, plus the gas phase free energy, that should be equal for all hemes, due to reasons already explained. Therefore, the free energy difference G from Fig. 2 a can be calculated entirely based on individual pseudoreduction processes performed for the hemes in question, in their protein environment. Consequently, the difference between redox potentials of different hemes can also be calculated.

The free energy associated with the individual processes described above can be calculated by using molecular dynamics techniques associated with free energy perturbation methods (Zwanzig, 1954) using the coupling parameter () approach (Mruzik et al., 1976), that drives the hamiltonean of the system from the oxidized to the reduced form (by changing the charges). "Artificial" thermodynamic integration (King, 1993) is used to compute the free energy difference between the oxidized and reduced state (Eq. 4, where it is considered that there are no changes in kinetic energy).

(4)

Improved method for free energy calculations with large electrostatic changes

Changes between entities that are quite different in terms of charges pose specific problems (Åqvist, 1990; Auffinger and Beveridge, 1995; Chipot et al., 1994; Straatsma and Berendsen, 1988). This is related to the long-range character of the electrostatic interactions, which extend beyond usual cutoff values and are therefore a source of error in this type of free energy calculations. One approach used to circumvent this problem is to use a Born-type correction, but this is only good when the medium beyond the cutoff can be approximated by a continuum, which is not true in the present application. In addition, this correction is sometimes inaccurate when short cutoff radii simulations are compensated by the use of the correction (Straatsma and Berendsen, 1988). Ewald summation is another method used, but is more suited for periodic systems. Besides these, there are a number of other methods to treat this problem (Smith and van Gunsteren, 1993). However, most of them are difficult to implement in current simulators and/or are not of universal applicability to all systems. We should refer that very recently, after the submission of this paper, an interesting approach to handling free energy calculations involving net charge changes appeared in the literature (Simonson et al., 1997). This method partitions free energy calculation into a process performed on a local system treated microscopically with free energy calculations, and into two other processes performed in the global system treated using continuum electrostatic calculations. In this way, the problem of the long-range electrostatic contributions for the free energy is transferred to the continuum electrostatic steps, which are less computationally demanding and can efficiently model the solvent (and protein) long-range dielectric response.

An obvious approach to diminish the errors associated with the cutoffs is to use a quite large cutoff. In the calculations of relative redox potential of azurin (Mark and van Gunsteren, 1994), an 18-Å cutoff was used, but this substantially increases the computer time required. In addition, the use of a large cutoff raises questions due to an altered behavior of the molecular system, especially the solvent, which can display different mobility and structure. In fact, the SPC water model (Hermans et al., 1984) was parametrized to reproduce water properties when a cutoff of ~10 Å is used. Given these problems, we propose instead the use of a method that, although not suited for correction of long-range electrostatic interactions in general, can correct the effect of short cutoff truncations in the calculation of free energies. This method is based on the usually called twin range method, but uses two, instead of one, long-range cutoff. The method is illustrated in Fig. 3.

View larger version (145K): [in this window] [in a new window]   FIGURE 3   The double long-range cutoff scheme for electrostatic perturbations. Within zone A all interactions are calculated every MD step. Within zones B and B', interactions are updated every n steps. For the majority of atoms the interactions are only calculated up to the zone B. However, for perturbed atoms, such calculation is extended over B'.

As in the normal twin range method, one considers one cutoff sphere (R1) where interactions are updated every integration step. Then, between R1 and R2, the interactions are calculated once every n steps (10 for instance) and in between, the force evaluated in this step is used. This allows for significant savings in computer time. Beyond R2, the interaction is not taken into account. This is a significant source of error in free energy calculations involving charged species: beyond the R2 cutoff (usually something between 8 and 14 Å) the interaction is still non-negligible. To improve this situation, another cutoff, R2', is introduced here, but this cutoff is only applied to interactions where at least one of the charge groups is part of the perturbation. In other words, charged groups involved in the perturbation see up until R2' and the other groups see only up until R2. Obviously, normal groups see the perturbed groups if they are at a distance R2', the electrostatic force always being applied in both directions. Since the groups involved in the perturbation are a minority, the use of this long cutoff does not substantially increase the computer time required to perform the simulation, but the results obtained can be quite different and improved in terms of electrostatic free energy. The calculations described in this work were carried out with R1 = 8 Å, R2 = 10 Å, and R2' = 20 Å. Note that the use of 10 Å for the maximum cutoff, as is used in a general case, would result in most of the charged groups capable of influencing the redox potential actually being out of the sphere of interactions, due to the large size of the heme group.

Besides the general problem described above concerning the electrostatic cutoff there is an associated problem, which is the existence of "noise" in the calculation of V/, the quantity necessary for the calculation of the free energy. This is illustrated in Fig. 4.

View larger version (17K): [in this window] [in a new window]   FIGURE 4   Reduction of heme IV. (a) Sampling at  = 0.0 started from an equilibrated simulation (see below) with cutoffs R1 = 8 Å and R2 = R2' = 10 Å. (b) Similar to (a) but with cutoffs R1 = 8 Å, R2 = 10 Å, and R2' = 20 Å.

It is obvious that the value of the derivative in Fig. 4 a follows four distinct distributions that are mixed along the dynamics. These distributions correspond to the entry and exit of charge groups from the cutoff sphere of the charge group involved in the reduction of the heme (iron and axial nitrogens). This clearly shows that there is a large discontinuity at 10 Å, and that beyond this value the interactions are far from being close to zero. Such a situation is disastrous for this type of calculation, since it is difficult to establish the necessary ensemble average of the quantity, and worse still, the obtained value may have no physical meaning. Contrary to this, when the special long cutoff for perturbed charge groups (R2' = 20 Å) is used, the result is completely different, as Fig. 4 b shows. The discontinuities observed before disappeared and there is only one distribution for the derivative. Therefore, it is much easier to calculate the ensemble average and this value has physical meaning. On the other hand, the fact that in the situation corresponding to the 20 Å cutoff there are no distinguishable distributions in the V/ quantity means that the effect of the groups beyond this cutoff must be small, given that these can enter and exit the cutoff sphere without noticeable effect.

Both from the theoretical point of view as well as in practice, the introduction of this improvement in the current free energy calculations is quite advantageous. Obviously, this should not be considered a general treatment of electrostatic interactions, but in the present case (and possibly for other types of free energy calculations) it is certainly an approximation that improves results without increasing the computer time required. In any case, it is superior to the use of a short cutoff. The method was implemented in GROMOS and we will use it in all the studies described hereafter, including equilibrations before free energy calculations.

	RESULTS AND DISCUSSION

Top Abstract Introduction Methods Results & Discussion Conclusions References

Standard molecular dynamics simulation in explicit solvent

Different initialization conditions can result in different trajectories of a molecular system, and it is therefore difficult to make conclusions based on a single simulation (Björkstén et al., 1994; Soares et al., 1995). Therefore, four 1-ns molecular dynamics simulations were performed starting from the same optimized structure, but with different sets of random velocities with a distribution at the temperature of 293 K. The four simulations were done at the same temperature and the results can be appreciated in Fig. 5, where the r.m.s. deviation from the x-ray structure is represented along the simulation.

View larger version (23K): [in this window] [in a new window]   FIGURE 5   Four molecular dynamics simulations of cytochrome c3 in water starting from different sets of random velocities. Open squares, run 1; open circles, run 2; open diamonds, run 3; filled diamonds, run 4.

The results show that the structure is not stable under these conditions, even if the simulation identified as run 4 shows a plateau in the r.m.s. deviation at the end of the run. Nevertheless, the r.m.s. is quite high, above 2.57 Å at 1000 ps. The large deformations observed are due to solvent penetration inside the protein. This is clear in runs 1-3. In these simulations, the -sheet from residues 8 to 21 and the loop from residues 52 to 62 "open" and expose the heme core, especially heme IV. In the case of run 4, heme IV moves downward (see the orientation of Fig. 1) and pulls part of the polypeptide chain of the protein downward. In this case, there seem to be lower deformations on the -sheet and the loop, which seems to be the reason why the r.m.s. is lower than in the other simulations.

Our main hypothesis to explain these deformations is that the heme potential is not good enough in its present form. This is not due to the modifications introduced in this work with the new potential (different interaction parameters and heme modifications), since simulations done with the old version of the force field [van Gunsteren and Berendsen, 1987; with Åqvist corrections on the water van der Waals parameters (Åqvist et al., 1993, 1994)] yielded quite similar results, with a r.m.s. deviation above 3 Å in long simulations (C. M. Soares, unpublished results). Given the instability observed in the simulations it was decided to perform some studies with a more gentle equilibration to find out whether this could be a possible cause. The protein-solvent interface can be quite "hot" at the beginning of the simulation and this can lead to a certain degree of deformation, as we have seen in other proteins (C. M. Soares, unpublished results).

Simulations done with a slow release (gentle equilibration) protocol

A more gentle procedure to start simulations that gave good results in other proteins (C. M. Soares, unpublished results) is to carry out a slow release of the protein in the already equilibrated solvent. This equilibration procedure consisted of performing molecular dynamics of the simulation system using harmonic positional restraints in selected parts of the system, but always keeping the solvent free. This consisted of 50 ps restraining all protein atoms with a force constant of 10⁴ kJ mol¹ Ų, followed by two 50 ps simulations where the same procedure was applied first to main-chain atoms and then to C atoms. Thereafter, the C atoms are released in three steps, of 10, 20, and 20 ps, where the force constant for restraining is 10³, 10², and 10¹ kJ mol¹ Ų, respectively. After that, everything was set free. Therefore, the initial 200 ps are only equilibration time. The results for two simulations with velocities generated with different random number generator seeds are shown in Fig. 6.

View larger version (18K): [in this window] [in a new window]   FIGURE 6   Two molecular dynamics simulations of cytochrome c3 in water using the slow release procedure. Squares, run 1; circles, run 2.

The results are quite similar to those obtained with the simple random velocities initialization. Therefore, the deformation effect does not seem to be due to the initial conditions, but to the force field itself. Some of the effects observed in the previous simulations are also observed here, namely the "opening" of the -sheet from residues 8 to 21 and the loop made of by residues 52 to 62, which results in a partial exposure of the heme core, mostly heme IV. The downward motion of heme IV is also observed, at least in run 1. In run 1 there is also one extra deformation in the N-terminal of the protein, a part contacting heme I (see Fig. 1), this contact being lost in this simulation. It is difficult to say whether this particular motion is an artifact or not, given that the N-terminal of a protein is by its nature (and by its specific topology here) usually more flexible than the rest of the protein. In our opinion, the main "destructive" effect of the simulation is the exposure of the heme core and the concomitant water penetration that starts to destroy the protein's fold.

To investigate the potential problems due to the heme parametrization it is advantageous to look at monoheme proteins. One such protein is cytochrome c₆ from Monoraphidium braunii (Frazão et al., 1995). This is a high-resolution structure of a 90 residue long protein with a c-type heme with one histidine and a methionine as the axial ligands of the iron. The reduced form was simulated (C. M. Soares, unpublished results) using the same conditions used in the cytochrome c₃ simulations. Under these conditions a simulation of 2.5 ns at 293 K was performed and the result is presented in Fig. 7.

View larger version (15K): [in this window] [in a new window]   FIGURE 7   A molecular dynamics simulation of cytochrome c6 in water using the slow release procedure.

The r.m.s. deviation in this simulation is substantially lower than in the cytochrome c₃ simulations, reaching here a value of ~2 Å. Most importantly, the structure seems to be stabilized after the first nanosecond and remains stable until the end of the run, a rather long period. Obviously, this simulation cannot constitute a definitive proof, since the protein is substantially different from cytochrome c₃. However, it is a good indication that monohemic proteins are more successfully described with the force field, probably due to the lower heme/protein ratio.

It is obvious that the present models are not sufficiently stable in long molecular dynamics simulations in order to provide an adequate treatment of equilibrium states, which is a prerequisite of free energy calculations. However, to look at the thermodynamics of redox changes in these systems, it may not be necessary to look at all dynamic aspects of this protein in solution (these, anyway, are not sufficiently sampled even in nanosecond simulations). From known examples in the literature, there do not seem to exist large conformational changes due to redox changes in cytochrome c (Berghuis and Brayer, 1992; Berghuis et al., 1994; Takano and Dickerson, 1981a,b) and in flavodoxin (Watt et al., 1991): different redox states seem to show remarkably similar structures. If the same is valid for cytochrome c₃, a partially restrained model, that can reduce the instability problems observed, may be adequate enough for our studies.

A partially restrained model to study the redox states of cytochrome c₃

A simulation model under these conditions was designed with the following characteristics: it is the same solvated system as the free simulations, but the C atoms of all residues, plus the iron atoms of the hemes, are position-restrained to their x-ray structure configuration, with a force constant of 10³ kJ mol¹ Ų. All protein and cofactor atoms (except iron) are allowed to move freely, as well as the solvent molecules. The model retains a considerable conformational flexibility, but it will not unfold as in the free simulations. Obviously, other common models could be used instead, such as the spherical extended wall region (van Gunsteren and Berendsen, 1990), but that would require a different model for each heme and, in addition, special forces would have to be applied to confine the solvent molecules around the region of interest. The model applied here uses standard periodic boundary conditions instead. The position restraint of the iron atoms was used to ensure the topology of the redox center, but its use may be questionable in these studies. Nevertheless, test simulations without the use of restraints on the iron atoms were carried out and they showed that the effect of the restraints is negligible in this case, both structurally and energetically.

As would be expected, the simulation of this model leads to a fast stabilization of the structure. After 300 ps of simulation the model is equilibrated and it can be safely used for equilibrium analysis. This will be the model representing the fully oxidized state of cytochrome c₃.

The simulation shows some interesting features concerning the structure and dynamics of the oxidized state. Hemes and axial histidines are found to be in similar conformations to those in the x-ray structure. The water distribution around the hemes presents particular characteristics, already observed in the crystallographic structure analysis of this cytochrome in particular (Matias et al., 1993) and cytochromes c₃ from other species and strains (Czjzek et al., 1994; Higuchi et al., 1984; Matias et al., 1996; Morais et al., 1995). The molecular dynamics simulation does show that some of the crystallographic waters are indeed associated with the protein in the simulated solution structure. In particular, this applies to the water molecules that make hydrogen bonds with the N1 atom of the axial histidines, when this atom is not making a hydrogen bond with the protein. For the period from 200 to 400 ps, Table 2 contains the hydrogen bonds maintained between the N1 atom of the axial histidines.

View this table: [in this window] [in a new window]   TABLE 2   Analysis of hydrogen bonds with N1 of axial histidines in the simulation of the oxidized state

It is obvious that the proton of the N1 of axial histidines is always involved in a hydrogen bond to an acceptor which can be either a peptide carbonil oxygen or a water. Even when this atom is somewhat exposed, as in His-106, which is bonded to heme IV, there is a hydrogen bond with a water molecule which becomes trapped at the surface of the protein. In general, these water molecules bonded to N1 are the most "fixed" water molecules found in the simulation, a characteristic that may be understood by their particular environment, which is generally shielded from bulk solvent and contains this particularly strong type of interaction.

The water molecules bonded to N1 are not the only water molecules fixed inside the protein during the simulation period. There are quite a number of water molecules trapped within cavities, although in certain cases there is exchange with the bulk water molecules. Some water molecules insert themselves between the heme plane and axial histidines (not making hydrogen bonds to N1). The majority is bonded to the main chain of the protein, but there are cases of hydrogen bonds with side chain atoms of the protein. Noteworthy are the hydrogen bonds between some water molecules and the propionate groups of some hemes, namely the propionates A and D of heme I and the propionate A of heme III. The presence of such water molecules may have structural and energetic importance.

The reduction of cytochrome c₃

The conformational and thermodynamic aspects of redox changes in proteins are still incompletely understood. The free energy perturbation approach already described may be a good approach to study the process of reduction, both in terms of free energy relationships as well as in terms of conformational effects. This is because the process of free energy calculation is itself a low perturbational way of introducing the redox change, and therefore safer when compared with an abrupt change. However, the path followed by the perturbational approach is not a physical path for this process, the only physically meaningful states being the initial oxidized and the final reduced states.

In the thermodynamic integration method it is necessary to sample the V/ function at selected values of  to be able to conveniently calculate the integral (the free energy). Therefore, it is necessary to know approximately the shape of the V/ function. This is done here with 100-ps slow growth simulations (where the change from oxidized to reduced states is coupled with the integration of the equations of motion), performing the reduction of each of the four hemes separately. The results of these studies reveal (data not shown) that the function V/ seems to be quite smooth and almost linear, which suggests that a uniform distribution of points in the thermodynamic integration might be the appropriate thing to do. This is confirmed by other longer slow growth simulations carried out for some of the hemes. To perform the thermodynamic integration calculations the following equally spaced values of  between 0 and 1 were used: 0.00, 0.25, 0.50, 0.75, 1.00. The simulations at each fixed value of  were carried out in a consecutive manner, i.e., the first was that with  = 0.00 (which was started from the structure at 600 ps of the equilibration simulation) and the last was that with  = 1.00. The last configuration of each sampling simulation was used as the starting structure of a slow growth simulation of 20 ps used to make a smooth transition between this  value and its next consecutive value. The samplings at fixed values of  were carried out over different amounts of time. For  = 0.00, 50 ps of sampling were used; for  = 0.25, 0.50, and 0.75, 100 ps were used; and for  = 1.00, 200 ps of sampling were performed. With this setup, each calculation of heme reduction took 550 ps total. Four of these calculations were carried, one for each heme.

Even with the use of a relatively long sampling time for each value of , there are some uncertainties in the convergence of the mean, given the large fluctuations displayed by the system. Actually, we have seen that other electrostatic perturbations with simpler systems like ions display the same type of behavior with large fluctuations (C. M. Soares, unpublished results). Despite these large uncertainties present in the convergence of the mean of V/, what is important is whether the difference between values of  is above the error associated with the sampling. This seems to be the case. The results of the thermodynamic integrations are presented in Fig. 8.

View larger version (13K): [in this window] [in a new window]   FIGURE 8   Averages of V/ along the perturbation. The standard deviation of the value is represented by error bars. (a) Heme I, (b) heme II, (c) heme III, (d) heme IV.

Despite the large fluctuations in the sampling of V/, which result in relatively large standard deviations, the V/ curve along  seems well defined with the chosen values of . The shape of the curves is different for each heme, hemes III and IV showing the largest deviation from a straight line. The free energy difference between initial and final states can be calculated by integration of the curves using the trapezoid rule. The free energy calculated in this way is presented in Table 3.

Reduction is thermodynamically more favorable for heme IV, in keeping with experimental results (Salgueiro et al., 1992; Turner et al., 1994, 1996). Therefore, the method has physical meaning, being able to predict ab initio the heme that is more prone to be reduced in the first step. The free energy difference for reduction of any heme at any stage of total reduction of the cytochrome can be derived from the experimental data (Salgueiro et al., 1992; Turner et al., 1994, 1996). Here we are interested in the free energy of reducing one heme at a time, keeping the others oxidized. The free energy depends on the pH, and this dependence was deconvoluted into a protonation effect (Turner et al., 1994, 1996). This protonation effect, which can represent one or two protons (Louro et al., 1996), is present and absent depending on the pH, with an apparent pK_a of 5.3 in the fully oxidized form, increasing to 7.4 in the fully reduced form. Therefore, from the experimental data one can derive two extreme situations, when the protonation effect is fully present and when it is not. Our calculated results can be compared with the two sets of results, and this is presented in Fig. 9.

View larger version (10K): [in this window] [in a new window]   FIGURE 9   Comparison of the calculated free energy with two experimental sets of experimental data of this system. (a) Comparison with results for the protonated form. (b) Comparison with results for the deprotonated form. The standard deviation of the value is represented by error bars.

It was explained above that free energy calculation methods are able to calculate relative free energies and not absolute free energies. In the previous plots, actual values of free energy are plotted, but the present analysis focuses on relative free energies. In fact, it is the correlation being analyzed, not the actual values. However, one should note the difference in the abscissa and ordinate axis scales. If the relative free energies were correct, the slope of the line fitting the points should be 1, since we are comparing the same quantity. Some possible reasons for the observed discrepancy will be presented below.

The results show correlation between the calculated free energy and the experimental results for the protonated form, while this correlation seems more or less absent when the results are compared with the experimental data for the deprotonated form. In Fig. 9 a, which corresponds to the comparison with the protonated form, despite the fact that heme II is lower than it should be (note also that the deviation is quite large), the results for the other hemes are almost in a straight line. Independently of the data used, reduction of heme IV is predicted the most favorable, a fact present in both forms. Therefore, despite the question of the energy scale, the method seems to have some physical meaning and predicts sufficiently distinct situations.

The lack of correlation with the experimental results of the deprotonated form, which should be the predominant form at pH 7, is surprising. This may have to do with specific protonation effects that are not taken into account here. This cytochrome is known to display a concerted proton and electron transfer (Louro et al., 1996), which may be associated with protonation/deprotonation dynamics of propionate groups of the heme (Soares et al., 1997), so this may be the cause of the observed dissimilarity. Two possibilities can be considered here: the first is to consider that protonation states of some residues are incorrectly assigned in the simulation. The second is to consider that the reduction process is actually more complex than simulated here, comprising not only reduction of individual hemes, but also simultaneous proton capture. The referred coupling can bring a new dimension to the thermodynamics of the process, and a correct treatment of these questions needs a simultaneous treatment of protonation states within the molecular dynamics simulation (Baptista et al., 1997).

The lack of correlation of the results for heme II is a puzzling observation. Again, specific protonation effects may be one possible explanation. However, conformational changes cannot be ruled out as an alternative, given that the constrained space used in the simulation may not be representative of the conformational space of the molecule in solution, at least in what concerns this heme.

Conformational effects of heme reduction

To assess conformational changes due to each reduction process, average structures obtained during the 200 ps prior and 200 ps after reduction were calculated for each heme. These average structures have information not only on the protein, but also on resident water molecules. To find these water molecules, water molecules closer than 3 Å from any protein atom at any given time were considered. From these, those whose average position displayed a variance <0.25 Ų were considered as resident waters. The data thus obtained allow the analysis of relevant changes due to the reduction of each heme. This is described in the following for each heme reduction.

Heme I

Heme II

View larger version (68K): [in this window] [in a new window]   FIGURE 10   Stereo pictures (Kraulis, 1991) showing the motion of propionate D from heme I as a result of reduction of heme II. (a) Before reduction. (b) After reduction.

Heme III

Heme IV

Contributions to the free energy of reduction

The reduction potential of a redox center can be controlled by a variety of factors (Bertini et al., 1997; Bertrand et al., 1995; Blackledge et al., 1996; Churg and Warshel, 1986; Gane et al., 1995; Gunner and Honig, 1991; Jensen et al., 1994; Lo et al., 1995; Mauk and Moore, 1997; Moore et al., 1986). A very important factor in controlling the redox potential of a redox center is the exposure to solvent. Solvent, in general, will have the effect of stabilizing the most heavily charged form of the redox center, therefore stabilizing the reduced state (1 for the oxidized state and 2 for the reduced state, including the charge on propionates). Preliminary results with a c-type heme model in water (a fully exposed heme) yielded a free energy of reduction of ~95 kJ/mol, while the same model in vacuum needs ~185 kJ/mol to become reduced. As expected, water stabilizes the reduced state when compared with the vacuum situation. This shows the large effect of solvent exposure in controlling the redox potential of a heme. However, solvent exposure is not the only factor that can play a substantial role in controlling the redox potential. Protein dipoles and ionizable groups are other factors to take into account. In general and in a simplified manner, the influence of negative charges will stabilize the oxidized form, while positive charges will tend to stabilize the reduced form. While the precise effect of specific charges on the redox potential of hemes would require further studies, the data obtained on heme reduction on cytochrome c₃ and on the heme model give us clues to the importance of this effect. The reduction of heme models yields positive free energies, which can be explained since the heme overall is electrostatically negative due to the presence of the propionates. However, the reduction of heme IV yielded ~135 KJ/mol, evidencing the effect of charged groups in controlling the redox potential, in this case by stabilizing the reduced state more than what simple solvent accessibility could do. In fact, around heme IV there is a patch of lysines (Matias et al., 1993) which can stabilize the reduced state by means of the positive potential that they create.

Critical analysis of the calculation model

The model presented here to study the effects of heme reduction and free energy calculation has physical meaning and allows the analysis of a number of interesting features. It is possible to look at dynamic effects and conformational changes due to heme reduction simultaneous with the calculation of the free energy. In addition, unlike the continuum electrostatic methods normally used to calculate redox potentials (Gunner and Honig, 1991), the methodology allows the analysis of problems where conformational changes are important for defining this redox potential. However, the method assumes a number of simplifications that may be responsible for its problems. Essentially, these are a result of the problematic treatment of electrostatic interactions in molecular dynamics simulation. Despite this treatment being greatly improved by the use of the double long-range cutoff described here and applied in the free energy calculations, there are obviously some missing or incorrectly treated factors. Manifestations of these are not only the inconsistency of the results for heme II, but most of all the different scale showed by the calculated free energies when compared with the experimental ones. First of all, despite the fact that the use of a long cutoff is needed for accounting for long-range electrostatic effects in the free energy, its use also represents an overestimation of the electrostatic interactions in the simulation: these are parametrized for using smaller cutoffs, indirectly trying to compensate for the long-range effects by increasing the short-range effects. Other sources of error in the model are related to dielectric effects, atomic polarizabilities, and specific ionic effects. In these simulations and since we are dealing with a microscopic model meant for simulations in solution, a relative dielectric constant of one was used. If this is a common choice in simulation, it is not necessarily the best one for the calculation of redox potentials. On one hand, the model assumes that the known dielectric effects inside the protein and in solution are going to be modelled microscopically by atomic dynamics, but it is known that simulations do not necessarily yield the usual values of the dielectric constant commonly used to treat the solution and the protein interior (Smith et al., 1993; Smith and van Gunsteren, 1994). On the other hand, atomic polarizabilities are not included and this may be an important factor in determining the redox potential. In continuum electrostatic models these effects are implicitly included in the dielectric constant, while in the PDLD model they are treated explicitly. Another important factor for the correct treatment of the electrostatic interactions determining the redox potential may be specific ionic effects. These are quite difficult to include in molecular dynamics simulations due to the limited time scale that is possible to simulate. Within this time scale, the ionic effects would be largely dependent on the initial positioning of the ions, and therefore our choice was not to include them at all, since this could be a bias. Nevertheless, their effect cannot be forgotten.

Finally, it should be mentioned that these free energy calculations dealing only with changes in charge may lose some of the features of the chemical process taking place. In fact, this simple methodology is not sensitive to fine details of heme electronic structure modulation by the protein. One relevant example may be the distortions observed in the hemes of c-type cytochromes (Hobbs and Shelnutt, 1995), a fact proposed to have importance in the modulation of the redox potential.

Improvements on this model for calculating redox potentials can be undertaken in the areas mentioned above. Another area that may be important is the development of new charge sets that can better characterize the reduction process. We feel that the heavily iron-centered redox change is an oversimplification. In general, improvements in the heme potential energy function are needed for a proper simulation of this system. In addition, the possibility of running longer simulations may help the proper parametrization of the force field on one hand, and on the other hand it would allow better sampling for the free energy calculations.

	CONCLUSIONS

Top Abstract Introduction Methods Results & Discussion Conclusions References

Long explicit solvent simulations of proteins with a large content of covalently bond heme groups, like the case of cytochrome c₃ reported here, are still difficult due to problems in parametrization. In fact, in the case of cytochrome c₃, the core of the protein is mainly formed by the four hemes and any problem in their parametrization and in the parametrization of the attachment to the polypeptide chain might lead to instability. That is what is observed with the current set of parameters used here, which means that improvements in this field are certainly needed.

Despite the problems reported with the free models, the restrained model used in this study proved to be adequate for the study of redox processes that are present in this cytochrome. This model, when combined with an improved treatment of electrostatic interactions in free energy calculations, allowed us to make ab initio semiquantitative considerations about the reduction thermodynamics of individual hemes. Conformational consequences of the reduction process can also be studied: we were able to see changes upon reduction which, despite being small, may be quite significant given the strong character of electrostatic interactions in low dielectric media such as protein cores. The strong and long-range character of the electrostatic interactions in the modulation of the redox properties of hemes is demonstrated by the need for a long-range cutoff correction. The use of this correction is not totally devoid of problems, but constitutes a semiquantitative model with negligible computational costs when compared with the standard method, allowing for improvements in the physical reality of the simulated system.